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ABSTRACT  OF  DISSERTATION 


A  Polynomial  Primal-Dual  Interior  Point  Method 
for  Convex  Programming  With  Quadratic  Constraints 

i>y 

Lee  James  Lehmkuhl 
Captain,  United  States  Air  Force 
Doctor  of  Science  in  Operations  Research 
The  George  Washington  University 

This  study  involves  the  solution  of  a  convex  nonlinear  programming  problem  using  a  primal-dual 
interior  point  algorithm.  The  problem  is  to  minimize  x  subject  to  0,  i=l  to  m,  where  each 

(j,(x)  is  a  concave  quadratic  fimction.  We  specify  certain  common  regularity  conditions  to  guarantee  the 
existence  of  a  solution  x*  and  a  vector  of  KKT  multipliers  u*  >  0  such  that  (x*,  u*)  solves  the 
associated  Wolfe  dual.  The  algorithm  is  motivated  by  the  following  optimality  conditions,  given  a  fixed 
par^uneter  p  >  0,  for  x  and  u  such  that  s,(x)  >  0  and  u^  >  0,  i  -  1,  •  •  •,  m: 

m 

VTj(x,ii)  ~  b  -  ® 

i=l 

Vi(^)  “  p  =  0,  i  -  1,  .  .  .,  m. 

The  algorithm  uses  Newton’s  method  to  approximately  solve  these  equations  for  a  decreasing  sequence 
{p*'}  where  p*|0  as  k  — »  +  oo.  A  step  size  procedure  maintains  feasibility  and  seeks  to  decrease  some 
merit  function. 

The  Newton  direction  obtained  has  several  interesting  features.  The  direction  in  the  primal 
variable  x  is  closely  related  to  the  Newton  direction  generated  by  a  Sequential  Unconstrained 
Minimization  Technique  employing  the  logarithmic  barrier  function  (SUMT).  In  fact,  the  primed 
direction  can  be  viewed  as  a  SUMT  Newton  direction  with  a  perturbed  barrier  function  Hessian  matrix. 
The  perturbation  depends  on  the  degree  to  which  u,j,(x)  p  0,  or  perturbed  complementary 
slackness  (PCS),  is  violated  for  each  i.  If  the  deviation  from  PCS  is  assumed  l.o  he  small,  the  primal 
direction  is  close  enough  to  the  SUMT  Newton  direction  that  it  may  be  used  as  the  search  direction  in 
a  SUMT  algorithm.  The  resulting  primal  variable  algorithm  retains  the  polynomial  computational 
complexity  shown  for  SUMT. 

The  Newton  direction  in  the  dual  variables  generated  by  the  primal-dual  algorithm  is  also  studied. 
A  modification  to  the  dual  step  size  and  direction  allows  the  progress  in  reducing  PCS  at  each  iteration 
to  be  calculated  and  therefore  controlled.  This  control  allows  the  deviation  from  PCS  to  be  kept  small 


throughout  the  primal-dual  algorithm’s  progress,  and  thus  the  polynomial  complexity  shown  for  the 
primal  variable  algorithm  discussed  above  follows  for  the  primal-dual  algorithm. 


(Completed  in  1993;  132  pages.) 
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This  study  involves  the  solution  of  a  convex  nonlinear  programming  problem  using  a  primal-dual 
interior  point  algorithm.  The  problem  is  to  minimize  b  x  subject  to  ^,(x)  >  0,  i-T  to  m,  where  each 
j,-(x)  is  a  concave  quadratic  function.  We  specify  cert<dn  common  regularity  conditions  to  gu=?rantee  the 
existence  of  a  solution  x*  and  a  vector  of  KKT  multipliers  u*  >  0  such  that  (x*,  u*)  solves  the 
associated  Wolfe  dual.  The  algorithm  is  motivated  by  the  following  optimality  conditions,  given  a  fixed 
parameter  p  >  0,  for  x  and  u  such  that  g-(x)  >  0  and  u,  >  0,  i  =.  1,  .  .  .,  m. 

m 

VL(x,u)  =  b  ~  j^u,Vj-(x)  -  0 
1=1 

UjJ,(x)  -  p  =  0,  i  =  1,  .  .  .,  m. 

The  algorithm  uses  Newton’s  method  to  approximately  solve  these  equations  for  a  decreeing  sequence 
{fn  }  where  p  |0  as  k  — ♦  -t-  oo.  A  step  size  procedure  maintains  feasibility  and  seeks  to  decrease  some 
merit  function. 

The  Newton  direction  obtained  has  several  interesting  features.  The  direction  in  the  primal 
variable  x  is  closely  related  to  the  Newton  direction  generated  by  a  Sequential  Unconstrained 
Minimization  Technique  employing  the  logarithmic  barrier  function  (SUMT).  In  fact,  the  primal 
direction  can  be  viewed  as  a  SUMT  Newton  direction  with  a  perturbed  barrier  function  Hessian  matrix. 
The  perturbation  depends  on  the  degree  to  which  u-jj(x)  -  p  —  0,  or  perturbed  complementary 
slackness  (PCS),  is  violated  for  each  i.  If  the  deviation  from  PCS  is  assumed  to  be  small,  the  primeil 
direction  is  close  enough  to  the  SUMT  Newton  direction  that  it  may  be  used  as  the  search  direction  in 
a  SUMT  algorithm.  The  resulting  primed  variable  algorithm  retains  the  polynomial  computational 
complexity  shown  for  SUMT, 

The  Newton  direction  in  the  dual  variables  generated  by  the  primal-dual  algorithm  is  also  studied. 
A  modification  to  the  dual  step  size  and  direction  allows  the  progress  in  reducing  PCS  at  each  iteration 
to  be  calculated  and  therefore  controlled.  This  control  allows  the  deviation  from  PCS  to  be  kept  small 
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througliout  the  primal-dual  algorithm’s  progress,  and  thus  the  polynomial  complexity  shown  for  the 
primal  variable  algorithm  discussed  above  follows  for  the  primal-dual  algorithm. 
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This  dissertation  involves  the  solution  of  a  convex  nonlinear  programming 

problem  using  a  prirnal-dual  algorithm  developed  in  McCormick(  1991a  and  1991b). 

The  problem  is  to  minimize  6Fx  subject  to  5j(x)  >  0,  i=l  to  m,  where  each  ^^(x)  is 

a  concave  quadratic  function.  We  specify  certain  common  regularity  conditions  to 

guarantee  the  existence  of  a  solution  x*  and  a  vector  of  KKT  multipliers  u*  >  0 

suck  that  (x*,  u*)  solves  the  associated  Wolfe  dual:  maximize 

L(x.u)  =  h^x  —  Y.  subject  to  V_cL(x,u)=0  and  u  >  0.  The  algorithm  is 

1=1 

motivated  by  the  following  optimality  conditions  for  the  problem  stated  above, 
given  a  fixed  parameter  //  >  0: 

Vj.L(xm)  =  0 

-  //  =  0,  i=L...,  m. 


The  algorithm  uses  Newton’s  method  to  approximately  solve  these  equations  for  a 
decreasing  sequence  {p*’}  where  as  k  — »  +  oo.  A  step  size  procedure  medntains 
feasibility  and  seeks  to  decrease  some  merit  function.  The  Newton  direction  is  the 


negative  inverse  of  the  Jacobian  matrix  of  the  above  equations  multiplied  by  the 
vector  of  the  equations. 

The  Newton  direction  obtained  has  several  interesting  features.  The  direction 
in  the  primal  variable  x  is  closely  related  to  the  Newton  direction  generated  by  a 
Sequential  Unconstrained  Minimization  Technique  employing  the  logarithmic 
barrier  function  (SUMT).  In  fact,  the  primal  direction  can  be  viewed  as  a  SUMT 
Newton  direction  with  a  perturbed  barrier  function  Hessian  matrix.  The 
perturbation,  that  is.  the  degree  to  which  the  primal  direction  differs  from  the 
SUMT  Newton  direction  depends  on  the  degree  to  which  \i^g-{x)  —  /x  =  0,  or 
Iierturbed  complementary  slackness  (PCS),  is  violated  for  each  i.  If  the  deviation 
from  PCS  is  assumed  to  be  small,  the  primal  direction  is  close  enough  to  the 
SUMT  Newton  direction  that  it  may  be  used  as  the  search  direction  in  a  SUMT 
algorithm.  The  resultinr  primal  variable  algorithm  retains  the  polynomial 
computational  complexity  shown  for  SUMT  by  den  Hertog,  Roos,  and 
Terlaky(1990).  The  Newton  direction  in  the  dual  variables  generated  by  the 
primal-dual  algorithm  is  also  studied.  This  direction  seeks  to  change  each  dual 
variable  to  adapt  to  changes  in  its  associated  constraint  5,(x)  as  x  changes  along 
the  primal  Newton  direction,  and  to  reduce  the  violation  of  PCS  at  the  current 
value  of  X,  A  modification  to  the  dual  step  size  and  direction  allows  the  progress  in 
reducing  PCS  at  each  iteration  to  be  calculated  and  therefore  controlled.  This 
control  allows  the  deviation  from  PCS  to  be  kept  small  throughout  the  primal-dual 
algorithm’s  progress,  and  thus  the  polynomial  complexity  shown  for  the  primal 
variable  algorithm  discussed  above  follows  for  the  primal-dual  algorithm.  Finally, 
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the  effect  of  the  modification  also  allows  the  use  of  the  results  in  Anstreicher(  1990) 
for  quadratic  programming  with  SUMT.  These  can  be  applied  to  show  another 
proof  of  polynomiality  for  linear  programming  with  the  modified  primal-dual 
algorithm. 
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CHAPTER  I:  INTRODUCTION 

I.l  Background 

Constrained  optimization  is  the  process  of  minimizing  (or  maximizing)  the 
value  of  a  function  of  one  or  more  variables  over  a  set  of  allowable  values  for  those 
variables.  The  set  of  allowable  values  is  known  as  the  feasible  region.  This  thesis 
presents  a  method,  or  algorithm,  of  constrained  optimization  for  a  particular  cl2iss 
of  problems.  Subsequent  chapters  contain  a  rigorous  definition  of  the  problem  to 
be  solved,  as  well  as  coverage  of  fundamental  theory  and  related  research;  but  for 
introductory  purposes  a  general,  non-rigorous  background  description  is  in  order. 
The  class  of  problems  considered  here  requires  the  minimization  of  a  convex 
quadratic  function  defined  over  E”,  known  as  the  objective  function.  The  feasible 
region  is  defined  by  one  or  more  concave  quadratic  functions,  known  as  constraints, 
in  the  following  manner;  the  feasible  region  is  the  set  of  those  values  of  x  G  E" 
which  cause  each  constraint  to  be  non-negative.  The  algorithm  moves  toward  the 
solution  along  a  series  of  points  in  the  interior  of  the  feasible  region,  that  is,  the 
subset  of  the  feasible  region  where  the  constraints  are  strictly  positive;  hence  the 
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term  “interior  point  method”.  The  problem  structure  provides  a  convex  feasible 
region,  which  together  with  the  convex  objective  function  results  in  a  convex 
quadratic  minimization  problem.  This  class  of  problems  is  of  significant  interest— 
for  instance,  a  constraint  on  the  maximum  Euclidian  distance  between  two  points 
may  be  expressed  with  a  concave  quadratic  constraint  like  that  mentioned  above. 
So  for  example,  the  minimization  of  the  distance  between  a  number  of  facilities 
subject  to  constraints  on  the  facilities  locations,  perhaps  within  a  certain  service 
area,  may  be  modeled  as  such  a  problem. 

There  are  several  important  classes  of  problems  which  conform  to  this 
structure.  Perhaps  the  most  common  and  well-known  is  the  linear  program,  in 
which  the  objective  function  is  linear  and  the  constraints  are  linear  affine.  Another 
is  the  convex  quadratic  program,  with  a  convex  quadratic  objective  function  and 
linear  affine  constraints.  Such  problems  occur  in  practice  and  require  solution.  For 
example,  linear  programming  may  be  used  to  plan  the  employment  of  a  limited 
arsenal  of  nuclear  weapons  to  maximize  the  damage  inflicted  on  a  foe,  or  to 
minimize  friendly  casualties  for  a  given  level  of  destruction.  Another  application  is 
the  allocation  of  limited  resources  (labor,  parts,  equipment,  and  so  on)  within  a 
factory  to  maximize  profit. 

Much  research  has  gone  into  developing  algorithms  to  solve  problems  of  this 
type.  Some  specialized  algorithms  apply  just  to  the  linear  program,  and  others  to 
the  convex  quadratic  program.  Initially,  the  main  goal  in  algorithm  development 
was  to  guarantee  convergence  to  an  optimal  value.  Then,  the  speed  of  convergence 
and  computational  considerations  came  into  focus.  Some  algorithms  which  perform 


well  in  practice  can  be  shown  in  a  worst  case  scenario  to  require  a  prohibitive 
amount  of  computation  effort,  while  others  which  possess  a  more  reasonable  upper 
bound  on  computational  effort  may  converge  so  slowly  in  practice  that  they  prove 
useless  for  solving  actual  problems.  One  can  thus  see  the  development  of  two 
measures  of  algorithmic  performance.  One  is  the  algorithm’s  speed,  in  terms  of 
computer  time,  in  solving  a  wide  variety  of  test  problems.  Another  is  the 
theoretical  bound  on  the  number  of  computations  required  in  the  worst  case  to 
solve  a  general  problem  of  a  given  size.  (The  measure  of  problem  size  will  be 
discussed  later.)  In  this  paper,  the  algorithm  presented  will  be  analyzed  for 
convergence  and  a  theoretical  bound  developed  on  computational  requirements. 

1.2  Organization  of  the  Dissertation 

The  dissertation  is  organized  into  eight  chapters.  This  introduction  is  Chapter 
I,  the  remainder  of  which  will  be  devoted  to  an  overview  of  the  contents  of  the 
subsequent  chapters. 

Chapter  II  covers  the  fundament£il  concepts  of  interior  point  methods  and 
reviews  some  current  research  in  the  area.  The  fundamental  concepts  portion  gives 
primary  emphasis  to  the  logarithmic  barrier  function  and  includes  a  simple 
example  in  E‘‘.  After  a  brief  discussion  of  unconstrained  minimization  with 
Newton’s  method,  the  chapter  introduces  the  concept  of  computational  complexity 
of  algorithms.  The  chapter  concludes  with  a  survey  of  recent  developments  in 
analyzing  the  computational  complexity  of  various  interior  point  methods. 
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Chapter  III  begins  with  the  rigorous  presentation  of  the  problem  to  be  solved 
and  the  Primal-Dual  Algorithm  which  is  the  core  of  this  research  effort.  Several 
lemmas  show  important  mathematical  properties  of  the  algorithm  which  will  be 
used  in  later  analysis. 

Chapter  IV  analyzes  the  search  direction  in  the  primal  variables  x  G  E” 
generated  by  the  Primal-Dual  Algorithm,  and  shows  it  to  be  closely  related  to  the 
search  direction  generated  by  an  interior  point  method  based  on  the  logarithmic 
barrier  function.  With  careful  bounding  of  the  difference  between  the  two 
directions,  a  primal  variable  algorithm  using  the  Primal-Dual  Algorithm  primal 
variable  search  direction  inherits  convergence  and  polynomial  complexity  results  of 
den  Hertog,  Roos.  and  Terlaky(1990)  for  the  logarithmic  barrier  function  method. 

Chapter  V  follows  the  primal  variable  analysis  of  Chapter  IV  with  similar 
analysis  of  the  dual  variable  search  direction.  An  attempt  is  made  to  show  a  type 
of  monotone  convergence  in  the  dual  variables  toward  stationarity,  but  examples 
show  that  such  monotonicity  may  not  occur. 

Chapter  VI  develops  a  modification  to  the  algorithm,  resulting  in  the  Modified 
Primal-Dual  Algorithm,  which  adjusts  the  direction  cind/or  magnitude  of 
movement  in  the  dual  variables  to  produce  the  monotonicity  sought  in  Chapter  V. 
Further  analysis  shows  that  the  Modified  Primal-Dual  Algorithm  inherits  the 
convergence  and  complexity  results  from  Chapter  IV. 

Chapter  VII  presents  an  alternative  convergence  and  complexity  result  for  the 
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Modified  Primal-Dual  Algorithm  based  on  the  work  of  Anstreicher.  By  applying 
the  Modified  Primal- L  ^al  Algorithm  to  a  linear  program  with  proper  initialization, 
we  show  convergence  and  polynomial  complexity  precisely  as  in  Anstreicher(1990). 

Chapter  VIII  discusses  techniques  for  meeting  the  initial  feasibility  and 
centering  conditions  required  by  the  algorithm. 

Chapter  IX  concludes  the  dissertation  with  a  summary  of  key  results  and 
potential  areas  of  future  research. 
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CHAPTER  II:  FUNDAMENTAL  CONCEPTS  AND  RELATED  RESEARCH 

II.  1  Fundamental  Concepts 

II. 1,1  Interior  Point  Methods 

This  section  will  present  some  basic  results  for  understanding  the  nature  and 
use  of  a  barrier  function  as  the  foundation  for  an  interior  point  algorithm.  The 
material  covered  is  not  intended  to  be  all-encompassing.  Readers  who  desire  more 
detail  are  referred  to  Fiacco  and  McCormick(1968,  1990),  from  which  much  of  this 
discussion  is  drawn. 

An  interior  point  method  as  employed  here  is  one  type  of  Sequential 
Unconstrained  Minimization  Technique  (SUMT).  As  the  name  implies,  the  idea 
behind  SUMT  is  solving  a  series  of  unconstrained  problems,  the  solutions  of  which 
converge  to  a  solution  of  a  constrained  minimization  problem.  There  are  a  variety 
of  Sequential  Unconstrained  Minimization  Techniques.  Some  allow  the 
unconstrained  problems  to  be  solved  on  the  boundary  or  outside  of  the  fea.sible 
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region  of  the  original  problem,  while  others  maintain  strict  feasibility.  The  later 
are  called  interior  point  methods.  Again,  a  variety  of  these  SUMT  interior  point 
methods  exist.  We  will  look  at  one  of  the  classic  types,  one  where  a  logarithmic 
barrier  function  is  used  to  enforce  feasibility,  and  a  gradual  relaxation  of  the  barrier 
leads  to  convergence  to  an  optimal  solution. 

To  be  more  precise,  let  us  define  the  general  constrained  minimization 
problem  as  follows: 

min  f(x) 

s.t.  ^^(x)  >  0,  i  G  {1.  2.  ...,  m} 

where  x  €  E”,  and  f  and  are  continuously  differentiable.  We  make  the  standard 
assumptions  of  the  existence  of  x®  such  that  >  0,  i  G  {1,  2.  ...,  m}  and  the 

existence  of  local  minima.  The  unconstrained  subproblem  is  defined  in  terms  of 
the  logarithmic  barrier  function  for  a  fixed  p  >  0: 

m 

min  $(x,^)  =  f(x)  - 

i=l 

s.t.  g-{x)  >  0,  i  G  {1,  2,  m}. 

The  logarithm  of  the  constraint  functions  forces  the  value  of  $(  •  ,/x)  to  become  very 
large  when  approaching  the  boundary  of  the  feasible  region,  so  minimization  of 
$(  •  will  move  away  from  the  boundary  (hence  the  term  “barrier  function'’). 

The  scalar  g.  is  known  as  the  barrier  parameter,  and  its  magnitude  determines  the 
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strength  of  this  barrier.  Hower-er.  for  any  //  >  0.  no  matter  how  small,  $(  •  ,/i)  will 
approach  +  oo  near  the  boundary  of  the  feasible  region.  In  practice  the  strict 
feasibility  constraints  are  essentially  ignored  when  minimizing  $(  •  .^).  since  the 
function’s  structure  itself  enforces  feasibility.  The  sequence  of  unconstrained 
minimization  subproblems  arises  from  minimizing  $(  •  for  a  strictly 
monotonicaay  decreasing  positive  sequence  It  is  shown  in  Fiacco  and 

McCormick  (Theorem  8)  that  the  solutions  to  these  subproblems  will  converge  to  a 
local  minimizer  of  the  constrained  problem.  Now  suppose  f(x)  and  —  ^,(x).  i  6  {1, 
2,  ....  m}  are  convex,  making  the  problem  convex.  Then  local  minimizers  are 
global  minimizers.  and  the  stibproblein  solutions  converge  to  a  global  minimizer. 

The  following  simple  example  shows  how  such  convergence  occurs  for  a  convex 
problem  in  E".  Consider  the  example  problem  from  McCormick(1983): 

min  X[  +  x.j 
s.t.  -  xj  +  Xo  >  0. 


For  /t  >  0  and  arbitrary,  the  logarithmic  barrier  function  for  this  problem  is 


<&(X,/X)  =  X^  +  X.,  -  -  xf  +  Xo)] 

and  its  gradient  and  Hessian  are 


V$(x,//) 


1  + 

2xi^ 

-  xf  +  X2 

1  -• 

•>  . 

-X1+X2 

8 


^^2 

0 


-  XT  +  X9 


[^J2  0]  + 


—  2x, 


X 


—  2xi 
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respectively.  On  the  interior  of  the  feasible  region,  the  Hessian  is  the  sum  of 
positive  multiples  of  two  rank  one  matrices.  Each  rank  one  matrix  is  the  product 
of  a  column  vector  with  its  transpose.  Therefore,  the  Hessian  is  positive  semi- 
definite.  and  since  the  two  vectors  are  linearly  independent,  the  Hessian  is  in  fact 
positive  definite.  So  =  0  within  the  interior  of  the  feasible  region  is  a 

sufficient  condition  for  a  strict  unconstrained  minimizer  of  $(  •  ,p).  Let  x(/r)  denote 
.such  a  minimizer.  This  condition  implies 


x(p)  = 


5  ’  1  +  ^ 


T 


For  instance,  if  =  1,  then  x(/i®)  =[“^1  Ip'-  Solving  a  sequence  of  these 
subproblems  for  {/i^}iO,  we  obtain  the  “central  trajectory”,  and  taking  the  limit  as 
yields  x*  =  “  5  ’  global  minimizer  to  the  original  constrained 

problem. 


Of  course,  most  problems  do  not  allow  an  analytical  solution  for  the  location  of 
the  central  trajectory.  The  solution  of  the  unconstrained  minimization  sub¬ 
problem  thus  becomes  a  significant  part  of  the  logarithmic  barrier  function  SUMT. 
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Much  research  has  gone  into  unconstrained  minimization  in  general,  resulting  in 
various  candidate  algorithms  for  minimizing  ^(  •  ,/z)  for  a  given  fj..  Steepest 
descent,  conjugate  gradient,  quasi-Newton,  and  Newton's  method  are  some  of 
possible  techniques.  Each  has  advantages  and  disadvantages,  and  a  detailed 
presentation  of  each  is  found  in  McCormick(1983).  However,  Newton's  method  has 
become  the  method  of  choice  for  the  unconstrained  sub-problem,  both  in  the 
developmental  days  of  SUMT  and  more  recently  as  the  speed  of  convergence  and 
computational  requirements  of  interior  point  algorithms  have  come  under  renewed 
study.  The  important  role  of  Newton's  method  in  underpinning  the  interior  point 
algorithms  studied  in  this  paper  calls  for  a  brief  discussion  of  the  application  of 
Newton's  method  in  unconstrained  minimization.  This  is  not  a  rigorous 
mathematical  development  such  as  is  found  in  McCormick(1983).  but  rather  the 
basic  ideas  needed  within  the  context  of  this  paper. 

The  initial  intent  of  Newton's  method  was  to  solve  a  system  of  n  equations  in 
n  variables.  The  Newton  search  direction  in  this  case  resulted  from  inverting  the 
.Jacobian  matrix  of  the  system  of  equations  and  multiplying  it  by  the  negative  of 
the  system  of  equations.  An  example  of  this  is  found  in  Chapter  III.  equations 
(3.4)  and  (3.5).  In  the  context  of  unconstrained  minimization,  the  gradient  and 
Hessian  of  the  function  to  be  minimized  replace  the  system  of  equations  and 
.Jacobian,  respectively.  Thus,  Newton's  method  is  used  to  find  a  point  where  the 
gradient  of  a  function  is  zero,  a  condition  known  as  stationanr\  "'.  Moriarity  is  a 
necessary  condition  for  a  local  minimizer.  If  the  function  is  -ri  unex  as  in 

the  previous  example,  this  becomes  a  sufficient  condition  for  n  ;!iinimizer. 

VVe  strengthen  the  continuity  assumption  to  require  f  and  —  i  £  {1.  2.  ...,m}  to 
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be  convex  and  twice  continuously  differentiable,  and  further  require  V“$(  •  ,/z)  to  be 
positive  definite  on  the  interior  of  the  feasible  region  (which  implies  strict 
convexity).  From  an  arbitrary  strictly  feasible  point  x^,  Newton's  method 
generates  the  next  point  in  the  search  for  stationarity  as  follows: 

x^  +  1  =  x'*-’  -(V2$(x^/i))' 

where  A  is  a  positive  step  size  scalar.  For  the  ‘"pure"  Newton’s  method,  A  is  equal 
to  unity;  otherwise,  the  value  of  A  is  determined  by  a  line  search  along  the  Newton 
direction  —  ( V“$(x^’,^))~  *V$(x^,/y)  using  $(  •  ,/i)  as  a  merit  function.  An 
important  way  of  motivating  Newton’s  method  may  be  seen  by  replacing  $(  •  .fi) 
with  a  quadratic  function  with  a  positive  definite  Hessian.  By  letting  A  =  1. 
Newton’s  method  finds  the  global  unconstrained  minimizer  in  one  step.  Thus, 
when  the  behavior  of  a  function  between  the  current  point  in  the  search  and  the 
unconstrained  minimizer  is  close  to  that  of  the  quadratic  approximation  of  the 
function  formed  at  the  current  point,  the  pure  Newton’s  method  will  perform  very 
well.  Indeed,  as  the  sequence  of  points  generated  by  successive  pure  Newton  moves 
nears  the  minimizer,  the  quadratic  approximation  becomes  better  and  better. 
Because  of  this  phenomena,  Newton’s  method  with  a  step  size  of  unity  is  said  to 
have  an  asymptotic  quadratic  rate  of  convergence  near  the  minimizer,  that  is,  the 
distance  from  the  current  point  to  the  minimizer  is  bounded  by  a  constant  times 
the  square  of  the  distance  from  the  previous  point  to  the  minimizer.  Another  form 
of  asymptotic  quadratic  convergence  for  Newtons  method  involves  the  magnitude 
of  the  Newton  step.  This  property,  which  plays  a  key  role  in  a  number  of 
complexity  results  for  interior  point  methods,  states  that  near  the  minimizer  the 
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magnitude  of  the  current  Newton  step  is  bounded  by  a  constant  times  the  square  of 
the  magnitude  of  the  previous  Newton  step.  However,  a  step  size  of  unity  may  not 
always  be  possible— such  a  step  length  at  times  may  lead  to  a  point  which  violates 
one  or  more  constraints  and  is  thus  not  allowable. 

Before  moving  beyond  Newton’s  method,  a  few  words  are  required  about  the 
linesearch  along  the  Newton  direction  when  A  =  1  is  not  used.  Bisection  or  golden 
section  search  techniques  are  among  those  used.  It  is  important  to  note  that  the 
strict  convexity  of  $(  •  ,/z)  on  the  interior  of  the  feasible  region  and  its  behavior 
near  the  boundary  applies  as  well  to  $(  •  ,/i)  evaluated  along  the  Newton  direction, 
so  a  strict  global  minimizer  with  respect  to  the  Newton  direction  exists.  It  will  be 
seen  later  that  for  theoretical  complexity  results  the  type  of  linesearch  is 
unimportant  as  long  as  the  computational  requirements  for  the  linesearch  are  of  the 
same  order  of  magnitude  as  the  computational  requirements  for  generating  the 
Newton  direction. 

II. 1.2  Computational  Complexity 

The  computational  complexity  of  an  algorithm  is  simply  the  theoretical  worst 
case  bound  on  the  number  of  arithmetic  computations  required  to  solve  a  given 
type  of  problem.  Such  a  bound  is  stated  as  a  function  of  certain  measures  of 
problem  size.  Three  such  measures  are  prevalent.  One  is  the  number  of  bits  of 
data  required  to  store  the  problem  in  a  computer,  and  is  denoted  as  L  for  length. 
The  following  method  for  computing  L  for  the  linear  programming  problem  {min 
c^x  I  Ax  >  b,  X  >  0}  is  found  in  Bazaara.  Jarvis,  and  Sherali(1990).  We  define 
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log(  • )  as  the  logarithm  to  the  base  2. 


L  =  fl  +  log(ii)  +  log(m)  +  E[1  +log(l  +  |c,|)] 

J 

+  E  E  [1  +  log(l  +  la,-  -I)]  +  E  [1  +  log(l  +  |b,|)]l 

I  j  i 

Another  measure  of  problem  size  is  the  dimension  in  Euclidian  space,  denoted  n; 
and  the  third  is  m.  the  number  of  constraints.  We  now  present  the  standard 
definition  of  the  function  0(  • )  used  to  express  computational  complexity,  using 
notation  from  Papadimitriou  and  Steiglitz(1982).  Let  f(n)  and  h(n)  denote 
functions  from  the  positive  integers  to  the  positive  real  numbers.  Then 

f(n)  =  0(h(n))  3  constant  k  >  0  such  that  for  large  enough  n, 

f(n)  <  kh(n). 

If  f(n)  =  0(h(n)),  we  say  “f(n)  is  of  order  h(n)”.  One  illustrative  example  of  the 
nature  of  0(  • )  is  that  for  any  constant  c,  cf{n)  =  0(f(n)).  An  expression  of 
computational  complexity  will  be  of  the  form 

N  =  0(h(n)) 

where  N  is  the  total  number  of  computations,  or  perhaps  iterations,  h  is  some 
function  from  the  positive  integers  to  the  positive  reals,  and  n  is  a  positive  integer 
or  a  vector  of  positive  integers  which  measures  the  problem  size.  For  a  given 
algorithm,  knowledge  of  such  a  function  h(  •  )  is  key  in  understanding  and 
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evaluating  the  complexity  of  the  algorithm,  because  h(  • )  provides  a  bound  on  the 
growth  rate  of  the  computational  requirements  as  the  size  of  the  problem  increases. 

As  an  illustration,  we  discuss  the  computational  complexity  of  the  simplex 
method  for  linear  programming  through  a  well-known  example  first  shown  in  Klee 
and  Minty(1972)  and  found  in  both  Bcizaara,  Jarvis  and  Sherali(1990)  and 
Papadimitriou  and  Steiglitz(1982).  The  gist  of  the  example  is  that  by  slightly 
perturbing  the  m  constraints  defining  a  m-dimensional  unit  hypercube  and 
minimizing  —  x^  over  that  perturbed  hypercube,  the  simplex  method  can  be  made 
to  visit  each  vertex  of  the  hypercube.  That  means  the  simplex  method  would 
require  2'”  iterations  to  solve  this  admittedly  worst-case  problem,  and  the 
computational  complexity  of  the  simplex  method  is  0(2*^).  Thus,  increasing  the 
dimension  of  the  problem  by  one  doubles  the  computational  effort.  This  is  an 
example  of  exponential  complexity,  from  the  exponential  function  K  >  1,  and 
indicates  that  the  computation  effort  is  multiplied  by  some  number  greater  than 
one  when  a  measure  of  problem  size  is  increased  by  a  constant,  usually  unity.  The 
implication  of  exponential  complexity  is  that  as  problem  size  becomes  large, 
computational  requirements  grow  so  fast  as  to  outstrip  the  capabilities  of  even  the 
fastest  computers.  It  is  important  to  state  clearly,  however,  that  exponential 
complexity  does  not  mean  an  algorithm  is  not  operationally  effective.  The  simplex 
method  is  a  good  example,  for  it  has  and  continues  to  solve  all  manner  of  real 
world  linear  programming  problems  quickly  and  efficiently.  The  exponential 
complexity  does  indicate  that  there  are  some  problems  which  may  take  a 
prohibitively  long  time  to  solve  by  the  simplex  method,  thus  prodding  researchers 
to  find  another  linear  programming  algorithm  without  this  drawback.  A  more 
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preferable  computational  complexity  would  be  expressed  by  a  polynomial  function 
of  problem  size,  for  instance  0(L''m®n^)  where  r,  s,  and  t  are  known  constants. 
Under  this  type  of  complexity,  the  increase  by  a  one  of  a  measure  of  problem  size 
results  in  the  addition  of  computational  work  of  lesser  order  (since,  for  example, 
0((n  +  1)^)  =  0(n^).)  Thus,  polynomial  complexity  avoids  the  “exponential  blow 
up'’  caused  by  multiplying  the  amount  of  work  by  some  K  >  1.  The  remainder  of 
this  chapter  will  exploie  recent  advances  in  polynomial  algorithms  for  linear 
programming  and  other  classes  of  constrained  optimization  problems. 


II. 2  Related  Research 

Recent  years  have  seen  a  plethora  of  polynomial  algorithms  for  various  classes 
of  constrained  optimization  problems.  This  section  will  present  an  overview  of 
several,  representing  some  different  approaches  to  solving  the  constrained 
optimization  problem.  Algorithms  for  linear  programming  are  the  most 
predominant,  for  at  least  two  reasons.  First,  there  exists  the  very  successful  and 
renowned  simplex  method,  albeit  with  exponential  complexity,  so  the  research  goal 
of  doing  “better  than  simplex”  is  tantalizing.  Secondly,  the  linear  affine 
constraints  allow  the  employment  of  linear  algebraic  techniques  (such  as 
projections  into  the  nullspace  defining  the  feasible  region)  not  suitable  for  more 
general  nonlinear  constraints.  This  attribute  applies  to  quadratic  programming  as 
well,  which,  though  not  as  ubiquitous  as  linear  programming,  is  also  heavily 
represented  in  the  literature. 
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These  algorithmic  overviews  are  intended  to  provide  the  reader  a  context  in 
which  to  view  the  complexity  results  presented  later  in  this  dissertation,  so  only  a 
concise  description  of  the  algorithm's  distinguishing  characteristics  is  warranted. 
Please  see  the  cited  references  for  complete  derivations.  When  an  algorithm  serves 
as  an  underpinning  for  this  dissertation's  results,  a  rigorous  development  of  that 
algorithm  is  included  at  the  point  w'hich  it  is  used. 

The  first  polynomial  algorithms  for  constrained  optimization  were  developed 
for  linear  programming.  The  first,  in  Khachian(1979)  and  the  second  in 
Karmarkar(1984)  are  both  covered  in  Bazaara.  Jarvis,  and  Sherali(1990),  from 
which  this  discussion  is  drawn.  Khachian's  algorithm  iteratively  generates  a  series 
of  shrinking  ellipsoids,  each  of  which  contains  the  set  of  points  solving  the  lineax 
programming  problem  defined  as  minimize  (c^x  |  Ax  >  b,  x  >  0}.  Khachian 
showed  that  the  algorithm  will  terminate  with  a  solution  after  a  number  of 
iterations  bounded  by  a  polynomial  of  problem  size,  and  if  that  bound  is  exceeded 
no  solution  exists.  The  computational  complexity  bound  of  Khachian's  algorithm 
is  0[(n  +  m)®L],  where  the  matrix  A  is  m  x  n  and  L  is  specified  measure  of  the 
input  length  of  the  problem.  Unfortunately,  the  actual  computational  performsince 
of  Khachian’s  algorithm  is  close  to  its  upper  bound,  while  the  simplex  method 
generally  achieves  a  complexity  bound  of  O(m^n)  in  practice.  Thus,  Khachian’s 
method  failed  in  application  when  compared  to  the  existing  simplex  method. 
However,  the  importance  of  Khachian’s  work  should  not  be  underestimated.  Like 
the  running  of  the  first  four  minute  mile  or  the  breaking  of  the  sound  barrier,  it 
showed  that  a  feat  of  previously  unknown  feasibility  could  be  done.  This  resulted 
in  a  resurgence  of  effort,  and  ultimately  the  “barrier”  of  polynomiality  was  broken 
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regularly. 

While  Khachian  broke  the  barrier  of  polynomiality,  Karmarkar  introduced  the 
first  polynomial  linear  programming  algorithm  which  also  competed  favorably  in 
practice  with  the  simplex  method.  Karmarkar’s  algorithm  addresses  a  particular 
form  of  linear  programming  problem:  minimize  {c^x  |  Ax  =  0,  e^x  =  1,  x  >  0} 
where  A  is  m  x  n,  A  is  of  rank  m,  n  >  2,  data  are  integer,  e  is  the  vector  of  ones, 
and  two  fairl}’  restrictive  assumptions  hold.  The  first  assumption  is  that  an  initial 
point  x^  =  (1/n,  1/n,  ...,  1/n)  is  feasible  and  the  second  is  that  the  optimal 
objective  function  value  is  zero.  (There  are  transformations  and  additional 
variables  which  allow  standard  linear  programming  problems  to  conform  to  these 
requirements.)  At  the  heart  of  Karmarkar’s  algorithm  is  the  following  projective 
transformation,  shown  here  at  the  start  of  the  (k  +  1)^^  iteration  with  a  feasible 

X*  >  0: 


y,  =  ^r^,i€{l,2,  ...,n}. 

The  algorithm  takes  a  step  in  this  transformed  space  to  minimize  the  (transformed) 
objective  function  over  the  intersection  of  the  simplex  constraint  e^y  =  1  and  a  ball 
centered  at  yg  =  (1/n,  1/n,  ...,  1/n)  with  a  radius  strictly  less  than  l/,Jn(n  —  1),  the 
distance  from  y®  to  a  facet  of  the  simplex  constraint.  This  minimization  is  easily 
accomplished  by  projecting  the  negative  gradient  of  the  transformed  objective 
function  onto  the  null  space  of  the  transformed  feasible  region  and  moving  along  it 
to  the  boundary  of  the  constraining  ball.  This  is  the  new  feasible  point,  ready  for 
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the  next  iteration.  After  O(nL)  iterations  and  computational  effort  of  O(n^’^L)  the 
algorithm  achieves  an  interior  point  solution  x  such  that  c^x  <  2  ~  at  which 
point  another  polynomial  algorithm  may  be  used  to  find  the  optimal  vertex.  (Such 
an  algorithm  will  be  described  later  in  this  section.)  This  bound  competes 
favorably  with  the  0(m  n)  practical  bound  of  the  simplex  method.  Of  particulax 
interest  to  logarithmic  barrier  function  research  is  the  potential  function  used  by 
Karmarkar  to  show  quantifiable  progress  at  each  iteration  toward  optimality.  That 
function. 


n/n(cx)  —  £  /n(x  ■). 

J-1  ^ 

is  a  type  of  logarithmic  barrier  function  for  linear  programming,  and  is  also  a 
component  of  the  method  of  analytic  centers  for  linear  programming.  This  aspect 
of  Karmarkar’s  algorithm  spawned  renewed  interest  in  interior  point  methods  that 
ultimately  lead  to  the  complexity  results  upon  which  this  research  is  based. 

Recent  analysis  by  Powell(1991)  of  Karmarkar’s  algorithm  with  the  potential 
function  shown  above  has  established  that  the  exponent  of  n  in  the  iteration  bound 
can  be  no  less  than  one.  This  is  shown  by  constructing  a  worst-case  linear  program 
in  with  n  inequality  constraints,  which  transforms  to  n  -f  3  variables  for  the 
equality  constrained  problem.  Karmarkar’s  algorithm  is  shown  to  terminate  in 
n/20  —  1  iterations  for  n  >  120  and  a  multiple  of  20. 

A  well-known  property  of  linear  programming  is  that  if  a  linear  program  has  a 
solution,  it  has  a  vertex  (or  corner-point)  solution.  When  applying  any  interior 


18 


point  method  to  linear  programming,  it  is  desirable  to  find  an  optimal  vertex  from 
the  nearby  interior  point  at  which  the  algorithm  terminated.  The  termination 
criterion  is  that  the  resulting  interior  point  x  satisfies  c^x  —  z*  <2~  z*  being  the 
optimal  value,  and  is  motivated  by  the  fact  that  for  any  non-optimal  vertex  x, 
c^x  —  z*  ^  (0  ,  2  “  Since  polynomiality  is  the  stated  goal  of  the  algorithms 
covered,  the  algorithm  for  finding  a  vertex  solution  should  also  be  polynomial  so  as 
not  to  negate  the  importance  of  that  property  for  the  linear  programming 
algorithm.  Such  an  algorithm  is  found  in  Bazaara,  Jarvis,  and  Sherali(1990)  and  is 
called  a  purification  algorithm.  It  is  based  on  concepts  first  set  forth  in  Charnes 
and  Kortanek(1965).  Gonzaga(1992)  also  discusses  the  basic  idea  of  a  purification 
algorithm;  “...  each  iteration  of  a  purification  algorithm  reduces  one  variable  to 
zero  along  a  descent  direction  for  cost,  doing  pivoting  like  the  simplex  method.  No 
more  than  n  iterations  are  needed,  with  0(n“)  computations  per  iteration 
(Gonzaga,  1992).”  Gonzaga’s  explanation  is  based  on  an  equality- constrained 
linear  program  in  E”,  and  it  is  clear  that  the  overall  complexity  of  this  technique  is 
of  O(n^).  In  practice,  getting  to  within  2  “  of  the  optimcd  value  may  be  very 
time-consuming  and  therefore  expensive.  A  practical  approach,  without  such  a 
stringent  accuracy  requirement  but  with  no  explicit  consideration  of  polynomiality, 
is  found  in  Mehrotra(1990). 

Another  established  linear  programming  interior  point  algorithm  is  known  as 
affine  scaling.  Gonzaga(1992)  introduces  affine  scaling  as  a  type  of  scaled-steepest 
descent  algorithm.  The  nomenclature  comes  from  the  scaling  of  the  problem  at 
each  iteration- the  scaled  space  is  achieved  by  dividing  x,-  by  the  current  point  x^, 
i  €  {1,  2,  ...,  n}.  This  maps  the  current  point  x^  to  a  vector  of  ones.  The  scaled 


objective  function  negative  gradient,  projected  onto  the  nullspace  of  the  scaled 
constraint  matrix,  is  the  descent  direction.  A  variety  of  step  length  procedures 
have  been  employed  along  the  affine  scaling  search  direction.  A  very  simple 
procedure  uses  a  constant  step  length  of  one.  Another  moves  a  fixed  percentage, 
say  95  percent,  of  the  distance  to  the  nearest  axis  to  maintain  strict  inferiority . 
Because  of  the  scaling  of  the  feasible  region,  affine  scaling  produces  an  ellipsoidal 
confidence  region  in  the  original  space,  as  opposed  to  the  spherical  confidence 
region  for  Cauchy's  steepest  descent  algorithm.  This  allows  better  progress  toward 
optimality  by  avoiding  the  "zig-zag"  prevalent  in  unsealed  steepest  descent.  See 
Gonzaga(1992)  for  a  more  thorough  coverage. 

Affine  scaling  for  linear  programming  also  exists  in  a  primal-dual  form.  An 
example  of  this  is  found  in  Huang  and  Kortanek(1991).  In  this  form,  the  search 
direction  derives  from  a  potential  function  which  incorporates  the  primal  variables 
(x)  and  the  dual  slacks  (s).  The  potential  function  is 

(p{x,s)  =  p  /n(x'^s)  - 

where  p  =  n  +  -fn.  The  authors  construct  primal  and  dual  search  directions  which 
reduce  the  potential  function  subject  to  primal  and  dual  feasibility  and  an 
ellipsoidal  constraint  on  the  scaled  magnitude  of  the  search  directions.  The 
resulting  algorithm  simultaneously  takes  steps  in  primal  and  dual  space  at  each 
iteration  to  solve  a  linear  program  with  computational  complexity  of  0(yfn  L). 

The  algorithms  discussed  so  far  have  been  for  linear  programming.  The 


algorithm  presented  in  this  paper  is  for  a  more  general  class  of  problem  with 
concave  quadratic  constraints  and  a  convex  quadratic  objective  function. 
Polynomial  algorithms  for  solving  problems  more  general  than  linear  programming 
problems  exist  and  are  well  documented,  so  the  remainder  of  this  chapter  will 
present  several  such  examples.  Of  course,  each  will  apply  to  the  more  restrictive 
linear  program  as  well. 

A  number  of  algorithms  for  solving  both  quadratic  programs  and  quadratically 
constrained  convex  programs  employ  the  method  of  analytic  centers.  General 
theory  concerning  analytic  centers  was  developed  in  Huard(1967).  and  an  example 
of  this  solution  technique  applied  to  quadratic  programming  is  found  in  Mehrotra 
and  Sun(1990).  Given  a  general  constrained  optimization  problem  in  E”, 

min  f(x) 

s.t.  x  e  IR  =  {^,(x)  >  0,  i  =  1,  2,  ...,  m}, 

for  some  z  G  let  the  region  IRj.  =  {x  G  E"  |  f(x)  <  z,  x  G  R}  be  bounded  with  a 
non-empty  interior.  The  analytic  center  of  is  the  point  in  R,  which  solves 

max  (z -f(x))'’n  ^,(x). 

i=l 

for  some  positive  integer  r.  An  equivalent  formulation  was  developed  in  Fiacco  and 
McCormick(1968,  1990)  and  is  related  to  the  logarithmic  barrier  function: 

m 

min  —  r/n(z  —  f(x))  —  /n(5r-(x)). 

i=l 
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Mehrotra  and  Sun(1990)  let  r  =  m.  a  common  practice.  With  this  definition  in 
mind,  their  method  of  analytic  centers  is  straightforward.  Let  f(x)  be  quadratic 
and  convex,  and  let  each  5j(x)  be  linear  affine.  We  are  given  z®  such  that  IR  q  is 
bounded  with  a  non-empty  interior,  and  x*^  a  close  approximation  to  the  analytic 
center  of  IR  q.  The  (k  +  1)^^  iteration  begins  with  z^  and  x*^  such  that  x^’  a  close 
approximation  to  the  analytic  center  of  IR  f..  The  algorithm  decrements  z^  to  z^  ^ 
and  finds  x^‘  ^ ,  a  close  approximation  to  the  analytic  center  of  IR  ^  j .  The  key 

behind  the  algorithm  of  Mehrotra  and  Sun  is  that  by  decrementing  z^  a  small 
amount,  the  algorithm  requires  only  one  partial  Newton  step  (a  step  size  less  than 
one)  to  arrive  at  an  approximate  analytic  center  of  IR  ^_j.  The  resulting 
computational  complexity  (including  the  Hessian  inversion  for  the  Newton  step)  is 
0(n^'/m  L). 

.Another  algorithm  using  analytic  centers  has  a  direct  bearing  on  the  algorithm 
developed  later  in  this  paper.  Jarre(1991)  presents  a  method  of  analytic  centers  to 
solve  a  convex  program  with  quadratic  constraints  and  a  quadratic  objective 
function,  precisely  the  same  problem  solved  by  this  dissertation's  primal-dual 
algorithm.  The  iteration  bound  is  0(>|m  |/n{e)))  for  a  tolerance  e,  also  the  same  as 
will  be  shown  for  our  primal-dual  algorithm.  The  interpretation  of  e  is  that  the 
algorithm  terminates  with  an  interior  point  x  such  that  the  difference  between  the 
objective  function  evaluated  at  x  and  at  the  true  constrained  minimizer  is  at  most 
e.  This  is  known  as  an  e-optimal  solution.  .Although  Jarre(1991)  uses  a  different 
formulation  for  the  definition  of  the  analytic  center  than  does  Mehrotra  and 
Sun(1990),  the  concept  is  the  same.  Given  a  small  decrement  of  z^,  the  analysis 
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shows  that  a  single  Newton  step  will  arrive  at  a  new  approximate  analytic  center. 
Jarre  also  briefly  discusses  the  problem  of  finding  the  initial  values  of  and  x*^. 
and  shows  that  by  adding  an  artificial  variable  and  constraint  the  same  polynomial 
method  of  analytic  centers  can  be  used  to  find  and  x^.  See  Chapter  V'lII  for  an 
explanation  and  adaptation  of  this  technique  for  meeting  the  initial  condition 
requirements  of  our  primal-dual  algorithm. 

The  emergence  of  the  user-defined  parameter  e  warrants  some  discussion. 
Because  it  is  user-defined,  e  is  not  related  to  problem  size.  .4  polynomial 
complexity  bound  containing  e,  such  as  that  of  Jarre(1991)  or  that  of  our  primal- 
dual  algorithm,  may  thus  seem  less  satisfying  than  those  for  linear  programming 
which  contain  only  problem  size  information,  namely  n,  m,  and  L.  The 
requirement  for  e  is  a  result  of  the  more  general  quadratic  objective  function  and 
constraints.  .4ssume  a  solution  exists.  With  a  nonlinear  convex  quadratic 
objective  function,  an  extreme  point  solution  is  not  guaranteed;  and  even  with  a 
linear  objective  function,  thought  an  extreme  point  solution  must  exist,  it  will  not 
necessarily  be  at  a  “vertex”.  Indeed,  the  concept  of  a  vertex  may  not  apply  in  the 
presence  of  nonlinear  quadratic  constraints.  For  example,  the  concave  quadratic 
constraint  g(x)  =  1  —  x^  —  y^  >  0  in  generates  a  feasible  region  consisting  of  the 
boundary  and  interior  of  the  unit  sphere,  which  is  compact  but  clearly  has  no 
vertices.  Thus,  there  is  no  way  for  some  measure  of  problem  size  to  indicate  a 
stopping  point  near  an  optimal  vertex  as  L  does  for  linear  programming.  Hence, 
the  use  of  a  tolerance  parameter  e.  One  must  take  care,  however,  that  the  manner 
in  which  e  appears  in  the  complexity  bound  does  not  compromise  polynomiality. 

As  a  reader  of  this  dissertation  pointed  out,  had  a  bound  of  0(\|m  2  ^  ^  )  resulted. 
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substitution  of  2  ~  for  e  as  is  done  for  linear  programming  would  result  in 
t'xponential  complexity  even  though  2  ^  ^  is  a  user-defined  scalar.  The  scale  of 

the  problem  also  may  complicate  the  selection  of  e.  For  an  arbitrary  e  >  0  and 
objective  function  f(x),  suppose  we  can  compute  T  such  that  |f(xj|  <  T  for  all 
feasible  x.  Then  with  the  simple  transformation  f(x)  =  (e/2M)f(x).  minimizing 
f(x)  over  the  feasible  region  results  in  every  feasible  x  being  an  e-optimal  solution. 
Jo  while  the  use  of  the  tolerance  e  is  necessary  for  the  result,  readers  should 
recognize  that  polynomiality  based  partly  on  e  does  differ  from  polynomiality  btised 
strictly  n  problem  size. 

The  research  wliich  provides  the  underpinnings  to  this  dissertation  concerns 
the  logarithmic  barrier  function  and  SUMT  for  solving  both  the  convex 
c  uadratically  constrained  problem  and  the  quadratic  program.  Significant  research 
ii,  the  recent  past  has  proven  polynomiality  for  this  solution  technique  for  both  of 
these  problem  classes.  Two  papers  in  particular  are  crucial  to  our  results--they  are 
mentioned  briefly  here,  and  covered  in  detail  later  in  subsequent  chapters  as  they 
apply,  den  Hertog,  Roos,  and  Terlaky(1990)  showed  that  the  classic  logarithmic 
barrier  method  solve.'  the  quadratically  constrained  problem  to  within  a  tolerance 
of  e  in  0(>lm  |ln(€)|)  iteravions.  The  convergence  and  complexity  analysis  of  den 
Hertog,  Roos,  and  Terlaky  provides  the  basis  for  the  initial  complexity  result  for 
our  primal-dual  algoritl’.m.  and  is  presented  in  detail  in  Chapter  IV. 
Anstrcicher(1990)  aoplied  the  logarithmic  barrier  method  to  the  quadratic 
program.  The  convergence  and  complexity  analysis  differs  substantially  from  den 
Hertog.  Roos,  and  Ter.aky(1990).  and  establishes  an  iteration  bound  of  0("Jrh  L). 
Furthermore,  Anstrcic  icr  shows  that  with  the  correct  initialization  and  selection  of 
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parameters,  the  classic  logarithmic  barrier  SUMT  as  implemented  by  Fiacco  and 
McCormick  in  the  1960s  has  polynomial  complexity  of  0('i|m  L  /n(L))  for  quadratic 
programming.  Anstreicher’s  results  are  presented  in  Chapter  VII,  where  they 
establish  that  the  primal-dual  algorithm  developed  in  this  dissertation  has  an  0(>lm 
L)  iteration  bound  when  applied  to  a  linear  programming  problem.  (This 
particular  bound  does  not  extend  to  quadratic  programming.) 

This  chapter  would  be  incomplete  without  citing  two  recent  and  extensive 
s  i  .ey  papers  on  interior  point  methods.  The  first,  Gonzaga(1992),  has  already 
been  cited  with  respect  to  an  affine  scaling  linear  programming  algorithm.  It 
contains  a  weal'h  of  information  on  primal,  dual,  and  primal-dual  interior  point 
methodo  for  linear  programming.  The  second  is  Wright (1992),  which  is  an 
excellent  primer  on  the  basic  theory  underlying  interior  point  methods  for  linear 
rnnd  nonlinear  programming,  with  special  emphasis  on  the  logarithmic  barrier 
function.  Wright  includes  a  polynomial  complexity  proof  for  a  linear  programing 
barrier  function  method  that  is  somewhat  similar  to  Anstreicher’s,  as  w'ell  as  a 
survey  of  primal  and  primal-dual  methods  and  linear  algebraic  issues  for  interior 
point  methodo. 

Fi  .ally,  in  almost  any  reference  concerning  barrier  functions  in  mathematical 
programming,  including  those  above,  one  will  find  citations  for  Fiacco  and 
McCormick(1968,  1990).  Their  classic  text  on  nonlinear  programming,  first 
published  in  1968  and  now  reissued,  remains  a  fundamental  source  for 
understanding  this  field  of  research.  Their  efforts  at  the  time  were  directed  toward 
nonlinear  programming  and  did  not  involve  computational  complexity,  but  did  lay 
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the  groundwork  for  much  of  what  has  followed  in  linear  and  nonlinear 
programming.  The  startling  and  satisfying  discovery  by  Anstreicher  of  the 
polynomiality  of  their  original  1960s'vintage  SUMT  for  quadratic  programming 
indicates  the  enduring  value  of  their  research. 
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CHAPTER  III:  THE  PRIMAL-DUAL  ALGORITHM 

III.l  Problem  Definition 


We  define  a  convex  problem  P: 

min 

s.t.  g-{x)  >  0,  i  €  I  =  {1,2,  .  .  .  ,m} 

X  e  E" 

where  6  is  a  vector  in  E”  and  each  5j(x)  is  a  concave  quadratic  function  for  i  €  I. 
Let  IR  =  (x  G  E"  I  5,(x)  >  (J,  i  6  1}  and  R*  =  (x  £  E”  |  y,(x)  >  0,  i  G  I}.  Let  uc 
assume  R  is  bounded  and  R*  is  non-empty;  and  let  us  define  the  following 
logarithmic  barrier  function  on  R*,  for  /x  €  E^  and  strictly  positive: 

^  -  2  ln(?,(x)). 
t=l 

The  gradient  and  Hessian  of  this  function  are  as  follows: 
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i=i 


Vg,(x)Vg.(x)-"  _  V‘^ff,(x) 

gii^^  9i(^)  _ 

Note  that  $(x,/x)  is  convex  in  iR*  for  fixed  /z  >  0  and  thus  has  a  positive  semi- 
definiie  Hessian.  This  structure  allows  for  the  following  lemmas; 

Lemma  III.l:  For  an  arbitrary  z  €  E”  and  nonzero,  and  x  £  R*. 

If  R  is  bounded,  then 


g^{x)V gi{x)'^z  =  0,  i  £  I 

(3.1) 

z'^V^5-(x)  z  =  0,  i  £  I 

(3.2) 

cannot  both  be  true. 

Proof:  Suppose  not.  Then  (3.1)  and  (3.2)  both  hold.  Let  i  £  I  be  arbitrary.  That 

(3.1)  holds  implies  Vg-{x)^2  =  0.  Now  let  y  =  x  +  Az,  A  £  E^.  Using  x  £  R*  and 

(3.2)  leads  to 

9iiy)  =  9iM  +  ^  z'^V^g-{x)  z  =  g^{x)  >  0. 

Since  i  was  arbitrary,  this  holds  for  i  £  I,  which  implies  y  £  R*.  But  jj  y  ||  -♦  +  oo  as 


i=l 
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A  -» -f-  oo,  so  IR’*'  is  unbounded.  This  contradiction  shows  the  supposition  to  be  false, 
proving  the  lemma.  □ 

Lemma  III. 2:  If  !R  is  bounded,  then  is  positive  definite  in  R*. 


Proof:  Let  z  €  E”  and  nonzero  be  given. 


z'^Vf$(x,/z)z 


z 


T 


9^i^)  I 


(3.3) 


_  z^'Vg-{x)Vg-{xfz\  ^  z'^(  -  V25^.(x))z'\ 

Since  5',(x)  >  0,  and  Vgj(x)V5j.(x)^  and  (  —  V^g^{x))  are  positive  semi-definite, 
i  G  I,  each  term  in  each  summation  is  non-negative.  Then  by  Lemma  III.l,  at  least 
one  term  must  be  strictly  positive,  proving  the  lemma.  □ 


The  specification  of  a  linear  objective  function  in  problem  P  is  not  restrictive.  Any 
constrained  optimization  problem  with  a  general  objective  function  r(x)  to  be 
minimized  can  be  transformed  into  a  constrained  optimization  problem  with  a 
linear  objective  function  by  defining  another  variable  z,  adding  the  constraint 
z  —  r(x)  >  0,  and  minimizing  z  subject  to  this  and  the  original  constraints.  Thus, 
given  a  problem 

min  f(x) 

s.t.  fi^j(x)  >  0,  i  e  I 
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where  each  is  a  concave  quadratic  and  f  is  a  non-linear  convex  quadratic 
function,  one  may  convert  it  to  the  form  of  problem  P  by  introducing  an  additional 
variable  z  G  and  using  the  transformation  above  to  achieve 

P^:  min  z 

s.t.  Pj(x)  >  0,  i  G  I 
z  -  f(x)  >  0 
X  G  E”,  z  G  E^ 

This  is  an  equivalent  problem  but  with  a  linear  objective  function  and  concave 
quadratic  constraints.  Note  the  implication  that  a  convex  quadratic  program  (QP) 
is  therefore  a  special  case  of  P. 

This  transformation  does  not  compromise  the  e-optimality  criterion  mentioned 
in  section  II. 2  and  in  Jarre(1991),  as  we  now  show.  Let  (x*,  z*)  solve  P^,  and  let 
(x*^,  z^)  be  an  e-optimal  solution  for  P..  By  construction,  x*  solves  the  original 
problem  and  z*  =  f(x*).  Since  (x^,  z^)  is  feasible  we  have  f(x^)  <  z^  and 
f(x^)  —  z^  <  0.  Then  using  the  definition  of  e-optimality, 

z^-z*<e 

z^  +  f(x^)-z^-z*<e 
f(x*)-z*<e 
f(x*)-f(x*)  <e 
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So  is  an  e-optimal  solution  to  the  original  problem. 


Since  P  is  a  convex  problem,  there  exist  first  order  sufficiency  conditions  for  a 
global  minimizer  for  P.  We  first  define  the  Lagrangian  for  P  using  the  notation 
u  =  (Uj,  U2,  ....  u,„): 

_  m 

L(x,u)  =  b^x  -  Eu-5-(x) 

Henceforth,  we  will  let  VL(x,u)  denote  Vj.L(x,u). 

First  Order  SnffiriencY  Conditions  for  P  (see  McCormick(1983)  and  Fiacco  and 
McCormick(1968.  1990)) 

3  X  €  R,  u  >0  such  that  VL(x  ,u  )  =  0  and  u^-  5-(x  )  =  0,  i  6  I 

X  is  a  global  minimizer  for  P 

We  also  have  the  associated  Wolfe  dual.  The  assumption  that  K*  is  non-empty 
satisfies  Slater's  condition  for  convex  programming,  so  (x  ,  u  )  solves  the  dual 
problem  D  (see  Fiacco  and  McCormick(1968,  1990)): 

D:  max  L(x,u)  =  b^x  —  X) 

i=l 

m 

s.t.  VL(x,u)  =  6  -  E  u,Vg',(x)  =  0,  u^  >  0,  i  G  I 
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Returning  to  the  logarithmic  barrier  function  as  defined  for  problem  P,  consider  the 
problem 


min  ^*(x,/i) 
s.t.  X  ^  R 

for  an  arbitrary  fi  >  0.  Let  x(/i)  be  the  minimizer  of  $(x,/z)  in  R*.  By  defining 
Ui(/u)  =  fj./ and  n(g.)  =  {ui(/i),  n^in),  .  .  .  ,  u^(/i)}  we  see  that  (x(//), 
u(^))  is  dual  feasible.  This  dual  feasibility  allows  the  formulation  of  a  lower  bound 
on  the  optimal  value  of  the  objective  function  b^x  (see  Fiacco  and 
McCormick(1968.  1990),  Theorem  22).  This  lower  bound  combined  with  the  upper 
bound  given  by  b^x(/x)  (since  x(/i)  is  feasible)  yields 

b’^x(/x)  -  <  b'^x  <  b'^x(p) 

j=l 

=>  b^x(/i)  —  mp  <  b^x  <  b^x(p). 

A  more  extended  development  of  duality  in  the  minimization  of  $(x,  p)  is  found  in 
Fiacco  and  McCormick(1968,  1990),  the  text  upon  which  much  of  the  preceding 
discussion  is  based.  An  important  lemma  relating  to  these  upper  and  lower  bounds 
when  minimizing  $(x,  p)  is  contained  in  den  Hertog,  Roos,  and  Terlaky(1990),  and 
is  presented  here  without  proof. 

Lemma  III. 3:  For  decreasing  values  of  p  >  0,  the  objective  function  b^x(p)  is 
strictly  monotonically  decreasing  and  the  dual  objective  function  b^x(p)  —  mp  is 
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monotonically  increasing. 


Proof:  Omitted.  See  den  Hertog,  Roos  and  Terlaky(1990),  where  the  authors 
prove  the  result  for  general  smooth  convex  programs.  The  first  part  was  first 
proved  in  Fiacco  and  McCormick(1968,  1990)  for  the  more  general  convex 
programming  problem.  □ 


III. 2  The  Primal-Dual  Algorithm 

This  research  concerns  convergence  and  computational  complexity  of  an 
algorithm  that,  applied  to  problem  P,  simultaneously  finds  the  minimizer  x  and 
the  associated  vector  of  dual  variables  (i.e.  KKT  multipliers)  u  .  The  following 
motivation  for  the  algorithm  is  from  McCormick(  1991b). 

Primal- Dual  Algorithm:  Solve  the  equations 

VL(x,u)  =  0  (3.4) 

^i9iM  =  i  e  I  (3.5) 

with  X  €  R*  and  Uj  >  0,  iGl,  for  a  positive  sequence  |0.  Satisfaction  of  (3.5) 
is  known  as  perturbed  complementary  slackness.  For  a  given  value  of  >  0,  let 
[x(/z^),  u(/z^)]  denote  a  point  satisfying  the  above  equations  with  x(/z^)  G  R*  and 
Uj(/z^)  >  0,  i  G  I.  Then  [x(/i^),  u(/x^)]  is  dual  feasible,  and 
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>  b‘'x  >  b^x(/')  -  2  u-(/')5-(x(/i'^’))  =  h'^x(n^)  -  m/'. 

i=l 

The  point  x(/x^)  is  also  a  minimizer  of  this  can  be  seen  by  substituting 

/ g^{x{fi^))  for  u-(/i*^)  in  VL(x(/z*^),  u(//^))  =  0,  which  results  in  V$(x(/x^),/i^)  = 
0.  $  inherits  the  convexity  of  problem  P,  so  stationarity  and  feasibility  are 
sufficient  for  a  global  minimizer. 


In  practice,  these  equations  would  not  be  solved  exactly,  but  approximately.  The 
basic  idea  of  the  algorithm  is  to  use  Newton’s  method  to  solve  (3.4)  and  (3.5)  with 
a  line  search  to  determine  step  size  and  retain  feasibility.  We  assume  that  at  the 
start  of  the  Newton  iteration  we  have  >  0  and  (x.u)  such  that  x  G  IR*  and 
Uj  >  0.  i  G  I.  Let  U  =  diag(Uj-),  G  =  diag(5'(x)),  and  Vg  be  the  nxm  matrix 
whose  i^^  column  is  Vg^{x).  For  simpler  notation,  let  g,  denote  The  Jacobian 
of  (3.4)  and  (3.5)  is 


J(x,u) 


V2L(x,u)  -Vp 
UVg'T  G 


Using  the  bordered  inverse  formula  in  McCormick(1983),  we  have 


M-1  M-'^VgG-^ 

-G'^UV^'^M-^  D 


where  D  =  G  ^  —  G  G  \  and  M  is  defined  as  follows: 


M  =  V2l(x,u)  +  Vg\JG-^Vg'^ 
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u. 


=  V-^(x,u)  +  E  V5.(x)^V.9,(x)-^ 
/;=i  yjV^J 


for  X  G  R*  and  u^-  >  0.  i  G  I.  Some  manipulation  reveals 


M  =  £u^ff-(x)j 
i=l 


5,(xr 


The  follo'wing  lemma  results  from  the  structure  of  M  and  its  similarity  to  the 
Hessian  of  $(x,/t). 


Lemma  III.4:  R  is  bounded  =>  M  is  positive  definite  in  R*.  for  u  -  >  0,  i  G  I. 


Proof:  Let  z  G  E”  and  nonzero  be  given. 


z^Mz  =  z^ 


m 


V^i(x)V^,(x)'^ 

V2#((x)' 

_ 1 

9ii^) 

Since  u^-  >  0,  g,(x)  >  0,  and  Vg^{x)S/g^{x)^  and  (  —  V^g^{x))  are  positive  semi- 
definite,  i  G  I,  each  term  in  each  summation  is  non-negative.  Then  by  Lemma 
III.l,  at  least  one  term  must  be  strictly  positive,  proving  the  lemma.  □ 


Note  that  for  positive  definite  M  and  //  >  0, 


(3.6) 
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/  =  1  L  ’ 

and  is  positive  definite. 

We  have  shown  that  for  bounded  R,  M  “  *  exists.  If  x  6  R*  then  G  ~  ^  exists.  So 
given  the  boundedness  of  R  and  interiority  of  x.  J  "  ^(x,u)  exists.  Note  that  the 
diagonal  nature  of  G  and  U  implies  that  the  effort  in  computing  J  ~  ^(x.u)  is  of  the 
same  order  as  in  computing  M  “  ^  that  is,  O(n^). 

This  explicit  expression  for  the  Jacobian  inverse  allows  the  formulation  of  the 
Newton  step  with  step  size  t.  Here,  e  is  a  vector  of  ones. 


M-‘ 


VL(x,u) 
UGe  —  ejx 


t 
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X 

_  VI  -- 1 

,2  V9,(x)/^j(x) 

+ 

L  2  =  1 

J 

u 

b^ 

TTl 

-fill  ^^i(x)/?i{x) 

+ 

G  ^efi  —  u 

_ 

2=1 

Now  let  p 


V5,(x)/^j(x) 


The  Newton  step  becomes 


x(t) 

u(t) 


P 

-UG-'V^'^p  +  G-^e^ 


We  next  explore  the  implications  of  the  Primal-Dual  Newton  direction  in  the 
primal  variables. 
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CHAPTER  IV:  PROPERTIES  OF  THE  PRIMAL-DUAL  NEWTON 
DIRECTION  IN  THE  PR.MAL  VARIABLES 


IV.  1  Dcvcloning  a  Primal  \'ariable  Algorithm 


Let  H  denote  We  have  noted  the  strong  connection  between  H  and  M/ n 

for  fixed  fi  >  0.  namely. 

Vff;(x)Vg,(x)'^  _ 

.9,(x)-  9ri^) 


VM  = 


;=1 


Indeed,  if  ’\-g-{x)  -  fc  =  0  for  all  i  G  I,  then  M/g  =  H. 


.\ovv  let  u.’  =  max 
i  G  I 


1  - 


b.<7,(x) 


Then  for  any  d  G  E'*. 


rn 

V9,(x)V,9-(x)'^ 

n  u,’2_^ 

;=1 

.'7,(x)'^ 

5,(x)  J 

d  < 


d'^'kl  d 
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<  d1 


2=1 


.9i(x)“ 


5/(x) 


which  implies 

(1  -  w)d'''H  d  <  d‘'i\I  d  <  (1  +  u;)d'^H  d,  (4.1) 

and  if  a;  <  1. 

■\|1  —  u;  lid  11^^  <  II d  11^  {'^•2) 

where  for  any  positive  definite  nxn  matrix  A,  ||dli^  =  (d^A  defines  the 
matrix  norm  of  d  with  respect  to  A.  We  will  now  adapt  the  results  of  den  Hertog, 
Roos,  and  Terlaky(1990)  to  develop  a  primal  variable  algorithm  using  the  primal 
variable  Newton  direction  p(x,u,/i)  =  —  (M/zi)-  ^V$(x,/i)  generated  by  the  Primal- 
Dual  Algorithm,  and  u)  sufficiently  small  to  obtain  convergence  and  polynomial 
complexity.  We  will  use  p  as  an  abbreviation  for  p(x,u,//)  when  clarity  does  not 
suffer  because  of  it. 

A  Primal  Variable  Algorithm:  This  algorithm  loosely  follows  the  SUMT  central 
path.  It  gets  close  to  the  central  path  by  doing  linesearches  seeking  to  reduce 
$(x,/i)  along  successive  Newton  directions.  When  used  in  this  manner  in  the 
course  of  a  line  .search,  ^’(x,/i)  is  called  a  merit  function.  The  Newton  direction  p  is 
the  one  generated  in  the  primal  variables  by  the  Primal-Dual  Algorithm.  We 
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assume 


ij  =  max 

i  6  I 


1 


< 


X 

62 


at  all  times.  This  assumption  is  not  blithely  made— it  is  quite  restrictive,  and  the 
reader  may  reasonably  wonder  whether  it  could  actually  be  satisfied  without 
sacrificing  polynomiality.  In  fact,  Chapter  VI  of  this  dissertation  presents  a 
method  for  doing  so. 


It  will  be  shown  that  this  bound  on  a;  keeps  M///  sufficiently  close  to  H  to 
allow  proofs  of  convergence  and  polynomial  complexity.  The  algorithm  proceeds 
by  doing  successive  linesearches  until  2u;>|l  +  u;  <  || p <t  =  g’Jl  —  uj  ,  a 
tolerance.  The  positive  lower  bound  on  this  tolerance  is  quite  small  because  of  the 
size  of  w,  so  if  it  is  approached  during  the  linesearch  the  upper  bound  will  be  easily 
met,  and  the  linesearch  may  be  terminated.  (For  u)  <  ■^,  the  upper  bound  is 
strictly  greater  than  the  lower  bound.)  When  this  tolerance  is  obtained,  w'e  reduce 
H  to  JI  =  (l  —  for  some  0  €  (0, 1).  The  algorithm  ultimately  achieves  an  e- 
optimal  solution  value  for  P.  A  point  x  is  said  to  be  e-optimal  if,  given  the  true 
constrained  minimum  z*,  6'^x  -  z*  <  e. 


Input 

Reduction  factor  0  G  (0, 1) 

Tolerance  t  =  i-^1  —  u;.  Assume  ta  =  max 

°  i  G  I 


1  - 


<  i  alwavs. 
~  d2 


T  is  an  accuracy  parameter. 
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G  R*  is  given,  as  is  /i^  <  1/T,  and  >  0.  i  G  I  such  that  j  p(x^.u^,/i*^)  j ^  <  r 


begin 


X  =  x'^,  u  = 
while  fi  >  T  12x0.  do 

begin  outer  iteration 

while  IIP  11^/^  <  r  do 

begin  inner  iteration 


A  =  arg^rnin^  J$(x,/z)  s.t.  (x  +  Ap)  G 


Pll.Vf/^  -  +  ‘^} 


end 


X  =  X  +  A  p 
end  inner  iteration 
H  =  {\-d)n 
end  outer  iteration 


IV.2  Convergence  and  Complexity  of  the  Primal  Variable  Algorithm 

Let  qx{x  +  d,  fi)  =  +  V^'(x,^)^d  +  ^  d^H  d,  the  quadratic  approximation  of 

$(x  +  d,  fi)  at  X. 

Let  q3.(x  +  d,  ^)  =  ^*(x,//)  +  V$(x,/i)^d  +  |  d^^M  d 
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=  $(x,/i)  +  V^(x,fz)'^d  +  i  dT(H  -  E)  d 


where  E  = 


=  qj.(x  +  d.  //)  —  ^  d^E  d 

Vff,(x)Vff-(x)'^ 


V^gz(x) 

^z(x) 


From  the  relationship  between  H  and  and  the  definition  of  uj  we  have 

qj(x  4- d. /z)  —  I  ujd  d  <  qj.(x  +  d. /z)  <  q3.(x  +  d.  /z)  +  ^  u,’d^H  d 

The  proof  of  convergence  and  polynomial  complexity  for  the  primal  variable 
algorithm  is  developed  through  a  series  of  five  lemmas.  The  synopsis  below  is 
based  on  one  in  den  Hertog,  Roos,  and  Terlaky(1990),  and  may  aid  the  reader  in 
following  the  development  of  the  convergence  and  complexity  proof: 

(1)  Lemma  IV.  1  gives  an  upper  bound  on  the  absolute  difference  between 
$(x  +  d,/z)  and  qj;(x  +  d,/z). 

(2)  Lemma  IV. 2  shows  that  if  a  proximity  criterion  based  on  the  Primal-Dual 
Newton  step  in  the  primal  variables  is  met  at  a  point  x,  then  x  lies  close  to  the 
minimizer  of  $(  •  ,/z).  “Closeness”  is  measured  using  the  norm  with  respect  to  the 
Hessian  of  $(x,/:z). 

(3)  Lemma  IV. 3  establishes  that  when  the  proximity  criterion  from  Lemma 
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IV. 2  is  not  met,  one  obtains  a  guaranteed  reduction  in  #(  ■  ,/i)  when  doing  a  line 
search  along  the  Primal-Dual  Newton  step  in  the  primal  variables. 

(4)  Lemma  IV. 4  gives  an  upper  bound  on  the  difference  between  and 

where  x{fi)  is  the  minimizer  of  $(  ■  ,/i),  when  the  proximity  criterion 
from  Lemma  IV. 2  is  met. 

(5)  Lemma  IV'. 5  gives  a  bound  on  the  absolute  difference  between  the  values 
obtained  by  evaluating  the  objective  function  of  P  at  x  and  x(/i),  when  the 
proximity  criterion  from  Lemma  I\'.2  is  met. 

With  our  previous  assumption  of  u;  <  ^  throughout  all  iterations,  we  can 
prove  these  lemmas  analogous  to  those  in  den  Hertog,  Roos,  and  Terlaky(1990). 

The  proofs  are  similar  to  those  in  the  referenced  work,  adapted  to  deal  with  the  use 
of  M/fj,  instead  of  H  in  the  quadratic  approximation  of  $(x  -t-  d,  /j,). 

Lemma  IV.  1:  Let  cj  <  ,  ildi^^ <  -Jl  —  w  for  an  arbitrary  d  E  E”,  and  let 

X  E  IR*,  /i  >  0  be  given. 

Then  x  4-  d  E  R*  and 

|$(x  +  d.  /i)  -  qi(x -t- d,  j  < 


Proof:  lldll^y.^  <  ^|1  —  u;  =>  lldll^  <  1.  Using  the  Taylor  expansion  of  $(x./x),  we 
have 
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$(x  +  d.  /z)  =  q^(x  +  d,  /z)  +  x;  t ■ 

i=3 


where  t  •  is  the  term  of  the  Taylor  expansion.  In  den  Hertog,  Roos.  and 

oo 

Terlaky(1990)  it  is  proved  that  for  lld||^  <  1.  x  +  d  G  IR*  and  XI  fj  converges  such 


oo 

E 

z=3 


< 


IldlP;^ 

3(1  -II  dll/,) 


Thus,  '5(x  +  d.  p)  =  q2.(x  +  d,  ^)  +  i  d^E  d  +  X! 

i=3 

$(x  +  d,  //)  -  q^(x  +  d.  n)  <  i  d'^u.’H  d  +  £  tj 

i=3 


lemma.  □ 


|$(x4-d, /z)-qa;(x  +  d, /z)|  <  ^u^ldll^d- 


IldiP// 

3(l-||d||//)’ 


proving  the 


Lemma  IV. 2:  If  u;  <  ^  and  2uz^jl  +  u;  <  ||p||^^^  <  ,  then  ||x  —  x(/z)  ||//  <  | 

IIP  Hi/ 


Proof:  llPlIyv//^  <  1^1  -  w  =>  l|Pl|^<|-  Let  h  be  arbitrary  such  that  11  h  11^  =  I 
llPli//- 


lp  +  li||ff  <  llPllw +lihlii/  -  fllPlI//  T 


(4.3) 
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Consider  values  of  $  evaluated  on  the  ellipsoid  {x  +  p  +  h  |  llhllj^  =  |  ||p||^}.  From 
Lemma  IV.  1, 

^x  +  p  +  \u^i)  >  q^(x  +  p  +  h. /i)  -  ^  MIp  +  h||tf  -  3( 

Since  X  +  p  minimizes  q^.,  Vq_j(x  +  p,  /i)  =  0  => 

qj,(x  +  p  +  h,  /i)  =  qj,(x  +  p.  /i)  +  Vq2.(x  +  p,  /i)'^h  +  ih'^V“q^(x  +  p,  fi)h. 

=  q^(x  +  p. /i) +1  h'^(iM)  h  (4.5) 

Substituting  (4.5)  and  (4.3)  into  (4.4)  yields 

$(x  +  p  +  h,  fi)  >  qx(x  +  p,  fi)+  ^h'^(iM)  h  -  5  w^llPllitf --7T - ^ 

J(1  j) 

>  qx(x  +  PiM)+  5  (1 -‘*^)llhli^  -  IIPII^ -^IPII^ 

>  q3.(x  +  p,  ^)+  5  (1  -  ^*^)f  (iPlI/f  -  llPlItf -^iPlIff 

>  qj.(x  +  p, /r)+  (I  -yi*^)||P||ff  -^iPllir 
^  =»  9  —  34a;  >  0,  which  with  ||p||^  <  |  implies 

<&(x  +  p  +  h,  p)  >  qj,(x  +  p,  p)+  (9  -34u;)||p|(^  “^iPlIff 
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>  q^(x  +  p, /i)+  (if  -34u;)llp||5/ 


From  Lemma  IV.  1  we  also  have 

$(x  +  p, /i)  <  q^(x+p, /i)  + 

Since  2u,'4m^  <  Ijp  ||,^ ^’  <  5  \\P\\ff  ,  we  have 
'hlx  +  p.  p)  <  q^(x  +  p. /i)  +  illPifc  + 

3(1  gj 

<  qx(x  +  p.  fi)  +  gfiipiif/ 

and  so  tj  <  ^  H  If  ~  ^(^  +  P?m)  <  $(x  +  p  +  h.  p). 

Thus  the  center  of  the  ellipsoid  has  a  value  of  $  less  than  any  point  on  the 
boundary.  By  applying  the  strict  convexity  of  $  we  see  that  the  minimizer  x(/i) 
must  lie  in  the  interior  of  the  ellipsoid,  proving  the  lemma.  □ 

The  requirement  that  2a;-vjl  +  u!  <  ||p||^/ does  not  significantly  complicate  the 
algorithm  since  u  is  fairly  small.  The  continuity  of  the  gradient  and  Hessian  imply 
that  if  in  the  course  of  the  line  search  this  condition  fails  (||P||^/^^  is  too  small),  one 
can  “back  out”  along  the  Newton  direction  to  find  a  point  that  does  satisfy  the 
conditions  of  the  lemma. 
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Lemma  IV.3:  If  u;  <  ^  and  ((p  ||^ >  g'jl  —  u,’  ,  the  decrease  A  in  $  after  a 
linesearch  along  p  using  as  a  merit  function  is 


A 


Proof:  Let  A  be  a  step  length  such  that  ||Ap||^  <  1  and  0  <  A  <  1.  From  Lemma 
IV.  1  and  the  definition  of  qx(x  +  Ap,  /z), 


$(x  +  Ap,  p)  <  q^(x  +  Ap. //)  +  5  A^u.’llpii]/ + 


3(1  -  A||p||^) 


=  «^(x, /z)  +AV$(x)"p  +  ^AV^Mp  +  |a2z.-||p|||/  + 


3(1  ~  A||pi|^) 


=>  $(x, /z)  -  $(x  +  Ap, /z)  >  —  AV$(x)'^p  —  i  A^p^jgM  p 


-  p  11^  - 


3(1  -  A||p(|^) 


Now 

-  V$(x)'^p  =  V$(x)'^(iM)  -  *  V$(x) 

=  V$(x)T(iM)  -  '(iM)(iM)  -  Vt(x)  = 

so 

$(x, /z)-$(x  +  Ap,  >  A||pi|2jy^  - 

-  iA^aziinii^  __4eA_ 
2^H|P|I^  3(1-A||p1|^) 


(4.6) 
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Now  let  A  —  —  .  Then  with  some  manipulation  (4.8)  becomes 

^llPll// 

‘5(x, /i) -«5(x  +  Ap, /i)  >  -^IPII//  -  M  -  tA? 


We  also  have  ij  <  ^  and||pi|^  > 


which  together  imply 


$(x, /i) -$(x  + Ap, /z)  > 

Thus  the  lemma  is  proved.  □ 

Corollary  IV. 3:  !££.<;<  ^  and  llP||jv//^  ~  -$(x  +  p,  /z)  >0 


Proof:  From  (4.8)  and  A  =  1  we  have 
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$(x, //) -‘5(x  +  p, /^)  >  5  IIP  11^  -u;jjPll5^- 


IPII// 


3(1-||P||//) 


ilip|& 


>  (2 -  u;  - ^) Up  11'^  >  0.  □ 


The  next  lemma  may  be  of  general  interest,  since  it  allows  one  to  bound  the 
difference  between  the  current  \'alue  of  $(x,/i)  and  the  minimum  value  $(x(/i),  /z). 


Lemma  IV. 4:  If  ^  and  2u,•^jl  +  <  11P||,\//^  <  |  \jl  —  u,’  .  then 

$(x,  /i)  -$(x(p),  p)  <4)1p||5/. 


Proof:  Let  d  =  x(^)  —  x.  From  Lemma  IV. 1  and  the  definition  of  q^., 

^(x(/z), /z)  >  qj(x  +  d, /z)  -  ^  Jdll^- 

=  $(x, /z)  +V$(x, /z)'^d  +  id^^Md  - 
=  ^(x, /z)  -p^^Md  +  id^^Md  - 

From  the  Cauchy-Schwartz  inequality  we  obtain  p^^M  d  <  ||p||^  Thus 

$(x(/z), /z)  >  $(x, /z)  -||Pl!^/^  ~  5‘^ldlltf  ~  3^  -ildll^) 

^x,  fi)  <  iiPiiM/^lldll^/^  -  + 
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<  (  1  +  u;  )  II  P  II  (1 (  1  —  )  il  d  j|jf|  +  :ja-j|  d  ||]i^  + 


IldlP^ 

3(l-lldll^) 


(4.9) 


where  (4.9)  follows  from  (4.2).  From  the  definition  of  d  and  Lemma  IV. 2  we 
obtain 


$(x, /x) -$(x(/x), /i)  <  §(1 +  u;)|(p||5^  -  ^lldlf^  +  Jdll^ +^^y4|j^(4.10) 

Again  applying  Lemma  IV. 2.  the  fact  that  —  ^  +  u;  <  0,  and  ||  p  ||^  <  we  have 
(4.11):  (4.12)  results  from  1  +  u-'  <  |. 

<t(x,  ;.)-*( x(m).m)<  i(l+--)l|P|l5, 

<  |(1 +“')|PI&  +  IPife  (4.11) 


(4.12) 


This  proves  the  lemma.  □ 


Corollary  IV. 3  further  shows 


^(x  +  P,  /i)  -  $(x(/x),  n)  <  $(x,  /x)  -  «'(x(/x),  /x)  <  4  IIPII^.  (4.13) 
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8 


Lemma  IV.5:  If  a;  <  and  2u;Njl  +  u;  <  ||  p  ||  < 

[b^x  -  b  ^x(/z)  I  <  /i 


then 


Proof:  (See  den  Hertog,  Roos.  and  Terlaky(  1990).)  At  x(p).  V$(x(p),  p)  =  0. 
Taking  the  inner  product  of  V‘L(x(/i)-  /O  with  (x-x(/i))  we  have 


b  ‘x-b'‘x(/i)  ^  '  _  n 


Now  is  concave  ^  g-(x)  —  (ij(x(n))  <  V5j(x(/i))^(x  x(/i)) 


b  'x  -  b  ‘  x(//)  -^/,(x(p))  ^  ^ 

f/7(x(w)) 


(4.14) 


We  apply  Lemma  IV. 4  to  obtain 


$(x, /i)-$(x(M), f;M5i(x))-M5i(x(/i)) 

1=1 


h^x-h'^x{g)  ^  (  5,(x))  ^ 

<4||P||5/<T6  (4-15) 


Define  u-  =  5i(x)/gi(x(/i))  and  w  =  ^  ■  Then  from  (4.14)  and  (4.15) 

W  -  'Eivi-  1)  =  W  -(  E  Vj )  +  m  >  0 

2  =  1  \«  =  1  f 


(4.16) 


;=1 


(4.17) 


and  combining  (4.16)  and  (4.17)  we  obtain 


rn  ,  m 

E  i\~m  <  »'  <  jg  +  E  /«(''/)• 

i=:l  !=1 


(4.1S) 


The  gooinetric-arithnietie  mean  inequality  implies 


\!=1  /  !=1 


Let  e  = 


ml2  *’(•  ^  ihen  obtain 

i=l 


_  1  I  '**  1  » 

mu  —  m  <  tu  <  +  In  f]  v-  =  tr  + 


j=l 


*  ~  !6 


<E  +  in/n(u  ) 


*’  ^  -  m  -  16m  ^ 


(4.19) 


■^ln\v  )  —  u  )  =  -= —  i  <0  for  u  >1,  and  —  1  >  0  for  0  <  u  <  1.  so  /e(u  )  —  u 
is  monotonically  decreasing  for  u  >  1  and  moaotonically  increasing  for  0  <  "r  <1. 


Suppose  u  >  1+  2;^ 


+  +1<  j^  +  („(i  +  2^V(l+2^) 


+  1. 
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The  third  order  Taylor  expansion  of  /n^l  from  1  is 


ln(  1+0 


1^-1 _ 1 


+  - for  some  q  e 


2^^my  2>fm  Srn 


1.  1+. 


2\|m 


The  last  term  is  maximized  bv  o  =  1.  so 


L  +  („P)_7,  +1  <  1  +  1  +  - 


16m 


1  +.  1 


16m  2^^m  Srn  24m'^^~ 


16m  24m^/“ 


rn  >  1  =>  ]^^  +  Irriv  )  —  V  +  1  <  0  which  contradicts  (4.19). 

So  the  supposition  is  false,  and  r  <  1  +  2^=  .  Now  suppose  0  <  TF  <  1  —  2;—  . 

We  have  shown  that 

ln(v  )  —  T  is  monotonically  increasing  in  this  region.  By  similarly  manipulating  the 
third  order  Taylor 

expansion  of  In^l  —  2;^^  we  again  can  show  +  /n(v  )  —  T  +  1  <  0.  producing 
the  same 


contradiction.  Thus  1 


2Nfm 


<  V  < 


1  +  2^^  which  leads  to  the  following: 


f  1 _ L-')-  1  <  w 

I,  ^  2Vm  J  -  or 


=>  —  ;7>Jm  <  lu 
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(V 

m 


_ 1-+— 1- 


_L  <  ^ 

2^fm  -  2'm 


=>  It’ <  ^^^m  =>  lu’i  =  ^b^(x  —  x(/t))^  <  i>Jrn 

|b^(x  —  x(|/))  I  <  .  which  proves  the  leinma.  □ 


Convergence  and  Comnlcxitv 

Theorem  IV.  1:  The  algorithm  requires  at  most  K  =  outer  iterations  to 

arrive  at  a  solution  within  an  f  of  the  optimal  solution  value  for  P. 

Proof:  fL^  =  (I  ~  From  duality,  we  have  b^x(//’)  —  /t^’m.  With 

Lemma  I\L5.  this  gives  us 

b^x^  -  z*  =  b^x^  -  b^x(/i^)  +  b^x(/i^')  —  z* 

<  |b^x^  —  b^x(/x^’)|  +  b^x(/i^)  -  z* 

<  +  fi^m  =  + 

So  we  need  a  value  of  K  such  that  (1  —  +  m  j  <  e.  Taking  logarithms,  we 

obtain 


K  > 


-  /n(e)  +  /n(i%fm  +  m)+ln{^^) 
-/n(l  -0) 


54 


We  have  assumed  <  1/e,  and  know  B  <  —In(l  —  B)  for  0  <  0  <  1.  This  implies 


-2/n(e)  +  /n(2m)  ^  - /n(e)  + /n(ivm  +  rn)  + /n(/) 
B  - 


Therefore  K  >  ^  —  2/n(t)  +  In{m)  +  ln[2))  satisfies  this  inequality,  proving 


K  =  O  i/n(f )  .  □ 


Theorem  IV .2:  The  total  number  of  inner  iterations,  say  N.  during  an  arbitrary 
outer  iteration  satisfies 


—  1  +  d-  ^lIn)  ,  where  6  —  is  from  Lemma  IV. 3. 


Proof:  Let  p  be  the  value  of  the  parameter  for  the  current  iteration,  and  p  the 
value  for  the  previous  iteration.  Recall /I  =  —  B)fi.  We  denote  the  sequence  of 

points  generated  by  line  searches  during  this  iteration  as  x^.  x“,  ...  x"'.  x^  is 
the  value  of  x  at  the  beginning  of  this  outer  iteration.  From  Lemma  IV. 3, 


$(x'^,  p  )  <  $(x^  p  )  -  N6. 


Also,  $(x,  -  .£  ln(p-(x))  =  ^  -  £  ln(p,(x))  + 


N(5  <  $(x0,  ^)-^>(x''^, 

r 
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<  ^>(x^  //)  -  ^l)  +  -  h^x^') 

Because  x*^  is  almost  centered  with  respect  to  /i,  we  have  $(x®.  pi)  —  $(x(/i),  p)  <  1. 
=>  $(x°.  p)-<I'(x^'.  /O  =$(x*^,  /i)-$(x(/i),  p)  +  $(x(p).  p)-^{x^\  p) 

<  \+^x(p).  p)-<^(x-^\  p) 

<  1 

=>  N<5  <  1  + -^(6'x^  —  6^x''^)  (4.20) 

Since  x*^  is  nearly  centered  for  p  and  x^  is  nearly  centered  for  Ji  ,  Lemma  IV. 5 
implies 


|b'^(x°-x(/z))| 

(4.21) 

|b'^(x^'  -x(7I  ))| 

(4.22) 

From  Lemma  III. 3,  0  <  b^x{p)  —  h^x{p  )  and  b^x{p)  —  mp  <  b^x^p  )  —  mp  . 

^  0  <  b'^x(p)  —  b^x(p  )  <  in[p  —  Jl).  (4.23) 

Combining  (4.21),  (4.22),  and  (4.23,  we  get 


56 


~b^x{^)  +  b^x{fi)  —  b^x{n  )  +  b^x{n  )  -  b'^x"'^ 

<  lfji'im  +  m{^ -JI ) +  in(l  -  (1  -  6*))  +  ^(1  -  ^)Vm) 

<  /i(’Jm  +  ^m) 

Substituting  in  (4.20)  vve  thus  obtain 

yS  <  \+-L(b'^x^~b^x-'^’)<l+j^{^m  +  em).[J 

Combining  Theorems  IV.  1  and  I\''.2.  the  total  number  of  inner  iterations  is 

|{H-p4^(Vin4-^m))0(i/n(f))  =  -h  Om))  0(ln(f)). 

Furthermore,  the  order  of  0  has  a  key  role  in  determining  the  bound  on  the  number 
of  iterations. 

^  ~  ^  iterations 

^  =  0(1)  =>  0^m|/n(e)|j  iterations 
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Since  each  iteration  requires  ()(n^)  computations  (the  Hessian  is  inverted),  the 
computational  complexity  is: 

I  ln{e)  ()  for  0  =  j 

O^n^m  j//i(g)j)  for  6*  =  0(1). 

Consider  the  special  case  where  all  constraints  are  linear  affine.  (This  implies  P  is 
a  linear  program.)  Then  by  letting  e  =  2“  ^  the  purification  algorithm  for  finding 
an  optimal  vertex  described  in  Chapter  II  and  in  Gon2aga(1992)  may  be  used. 
Since  the  complexity  of  the  purification  algorithm  is  the  same  as  that  of  the 
Hessian  inversion,  the  overall  complexity  of  the  algorithm  is  not  altered. 
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CHAPTER  V;  THE  PRIMAL  DUAL  NEWTON  DIRECTION  IN  THE  DUAL 
VARIABLES 

V.l  Analysis  of  the  Dual  Variable  Newton  Direction 

In  the  previous  chapter  we  found  that  when  the  deviation  from  perturbed 
complementary  slackness  is  kept  small,  a  primal  variable  algorithm  using  the 
primal  portion  of  the  Primal-Dual  Newton  direction  will  converge  to  within  an  e  of 
the  solution  to  problem  P  in  polynomial  time.  It  is  therefore  logical  to  consider 
next  the  changes  in  the  dual  variables  induced  by  a  Primal-Dual  Newton  move. 
Indeed,  if  the  dual  variables  change  in  a  manner  which  maintains  the  small 
deviation  from  perturbed  complementary  slackness,  the  polynomality  of  the  Primal 
Variable  Algorithm  will  also  apply  to  the  Primal-Dual  Algorithm.  In  this  chapter 
we  investigate  the  nature  of  the  dual  variable  movement  and  how  it  affects 
perturbed  complementary  slackness.  It  will  be  shown  that  the  dual  variable 
Newton  direction  with  a  standard  step  size  procedure  does  not  guarantee 
maintenance  of  a  small  deviation  from  perturbed  complementary  slackness  in  all 
dual  variables.  However,  further  analysis  yields  a  modification  to  the  step  size 
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which  does  maintain  this  small  deviation. 


Recall  the  Primal-Dual  Newton  move: 


x(t) 

u(t) 


u 


-t- 


p 

-UG-'V^'^p  +  G-^e^l 


t 


u 


where  p  =  —  M 


m 


/=1 


=  -(^i) 


-  1 


V$(x,p). 


Let  i  G  I  be  arbitrary,  and  u-(t)  denote  the  new  value  of  the  associated  dual 
variable  given  step  size  t.  Then  we  have 


u-(t)  =  u.  + 


t 


^  -  “O'- 

The  first  term  of  this  expression  shows  how  Uj(t)  is  influenced  by  tp,  the  Newton 
step  in  the  primal 


variables.  ( 1 - r^VoTxl^tp 

9^)1 9ii^  +  tp). 


is  a  first  order  Taylor  series  approximation  of 


Sfj  is  concave  =>  l/^j  is  convex;  and  since  g',(x)  >0  =>  g^{x)lg-{x  -f  tp),  a  function 
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of  t,  is  convex. 


5,(x  +  tp)  -  r7-(x) 


So  multiplying  u^-  by  ^  ^tp  j  attempts,  using  first  order  information, 


to  account  for  the 


change  in  as  x  goes  to  x  +  tp.  keeping  Uj^j(x)  and  Uj(t)5^(x  +  tp)  approximately 
equal.  However,  the  inequality  above  shows  the  approximation  may  underestimate 
9ii^)/9ii>^  +  ^P)~  that  is. 

-^:|^^.9i(x)  '’tp)5,(x  +  tp)  <  u,^-(x). 

This  underestimation  may  be  observed  by  assuming  a  linear  affine  ^^(x). 

=>  g-{x  +  tp)  =  g'{x)  +  V^j(x)^tp.  Substituting  into  the  formula  for  nj{t)g-{x  +  tp), 

-  ^^^?i(x)’^tp)9.(x  +  tp) 

=  u.5,(x)(^l-^^V5-(x)'^tpj2^ 

=  Ui5i(x)  -  u-5-(x)|^^V^-(x)'^tp ^2 

Since  Uj  >  0  and  ^^^(x)  >  0,  the  last  term  is  always  non-positive,  and  (assuming 
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t  >  0)  equals  zero  only  if  Vg-{x)  =  0.  So  for  linear  affine  5,(x),  the 
underestimation  occurrs  whenever  S7g-{'x.)^-p  ^  0. 


More  will  be  said  later  about  how  this  underestimation  could  affect  the 
I)erformance  of  the  Primal-Dual  Algorithm,  but  for  now  let  us  move  on  to  the 
second  term  of 


To  study  the  role  of  the  term  ~  achieving  perturbed 

complementary  slackness,  we  let  p  =  0  =^>  g^{x  +  tp)  =  g,(x). 


+  tp)  -  n 


-/i 


=  (1  -t)u^-  + 


fit 

Fi(x) 


=  (1  -t)u,-5,(x)  -(1  -t)/r 


=  (1  -t)(Ujffi(x)-p) 


The  second  term  seeks  to  decrccise  |u,-5rj(x)  —  /i|  by  a  factor  of  (1  —  t);  t  =  1  reduces 
|uji9j(x)  —  /r|  to  zero.  But  this  is  predicated  on  the  first  term  adjusting  the  value  of 
Uj(t)  to  account  for  the  change  in  g^,  keeping  Uj(t)5fj(x  +  tp)  =  u-g^{x).  Thus  the 
terms  work  together;  the  second  term  attacks  the  current  deviation  from  perturbed 
complementary  slackness,  and  the  first  adjusts  the  value  of  Uj  such  that 
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Uj(t)5'j(x  +  tp)  =  n-g-{x)  so  as  not  to  confound  the  effect  of  the  second  term. 


V.2  The  Effect  of  the  Dual  Newton  Direction  on  Convergence 

To  study  how  the  dual  variable  Newton  step  performs  in  reducing  juj5j(x)  —  pj, 
we  use  the  value  of  the  step  length  t  introduced  in  the  proof  of  Lemma  IV. 3, 
t  =  l/(9l|Pl|;/)-  Besides  guaranteeing  a  reduction  in  ‘L(x.  p),  this  value  of  A  has 
other  valuable  properties: 

f  =  l/(9liP||//)  l|tp||//  -  5  =>  tp'^Htp  =  i.  (5.1) 

The  definition  of  H  implies 
tp  ^Htp  =  tp  ‘ 

All  the  matrices  in  this  summation  in  (5.2)  are  positive  semi-definite,  so  for  an 
arbitrary  i  €  I,  we  have 


gy5i(x)V5j(x)'’ 

t=l 


5,(x) 


2 


V“gi(x) 

9i(^) 


tp  =  ^. 


(5.2) 


5e(x) 


'^5*(x)’^tp 


(5.4) 
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(5.5) 


tp^ V"g.(x)  tp 

5,(x)  “ 


0  < 


rp 

2g,(x) 


< 


1 

162 


(5.6) 


0  < 


tp'^V“g-(x)  tp  j_ 
Qii^)  " 


(5.7) 


We  are  given  that  g,(x)  is  a  concave  quadratic  function.  Using  the  definition  of 
Uy(t)  and  a  2  °  Taylor  series  expansion  of  g-{x  +  tp), 


'h'(t)g,(x  +  tp) 


■^?,(x)^tp  + 


tp’^V“.9i(x)  tp\\ 

■29  M)  )) 


+  u,g^(x) 


Now  define  Wj  =  1 - ^  ^  u^g^(x)  —  fi  =  g,[—  u;^).  (Note  u)-  <  1  always.)  To 

study  the  reduction  in  |ujgj(x)  — /xj,  consider  two  cases: 


Case  ii  u-g,(x)  —  g  >  0.  This  implies  u;-  <  0.  Using  the  inequalities  in  (5.3) 
through  (5.7)  and  dropping  terms  with  the  appropriate  sign,  we  develope  the 
following  inequalities. 

Ui(t)g^(x  +  tpj -/X  >  (Uig-(x)-/x)(l-t(l+i))  +  u-g.(x)(--p^-^) 
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>  -  u.-)(l  -  -^t)  +  /i(l  -  a;-)(  - 

Iterations  only  occur  when  IIP II//  >  i  so  t  =  l/(9||p||// )  ^  t  <  §•  Thus. 
u,.(t)5.(x  +  tp)-p  >  +  + 

which  is  clearly  a  negative  lower  bound.  Now.  does  this  value  of  t  guarantee  a 
reduction  in  —  //? 

Ui(t)5i(x  +  tp)-p  <  M-s)(l-t(l-5-Tfe))  +  /^(l--'t)(Ti8) 

<  p(  -  + 

To  be  certain  of  a  reduction  in  u-g^{x)  —  p,  we  need  to  show 

p( -|||t  +  j^)  +  p(j^)  <  p(-u;-) 

^  ^i\  162'^  1458/  ^  1458 

If  0  <  —  i^i<  we  need  -  j||t  +  <  —  j||g,  or  t  >  0.048951.  So  if 

IIP  II//  >  2.27,  setting  t  =  l/(9||p||//)  may  not  give  a  reduction  in  u-g-[x)  —  p. 
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Case  2:  u^y^(x)  —  //  <  0.  This  implies  >  0.  Using  the  iiKMpialities  in  (5.3) 
through  (5.7)  and  dropping  terms  with  the  appropriate;  sign.  w('  (h'veiop  the 
following  inequalities. 


-  /'(  ~  “  irO  ^  ~ 

<  /d  ~-,d(i  -T^  +-nfe)  +  ^(ii) 

-Us  hefore.  w<'  know  t  <  |^.  Thus. 

s  sdnW) 


which  is  a  cpiite  small  positive'  upper  be)und  for  reasemabh'  values  of  //.  Xow.  does 
this  value;  of  t  guarante-e'  a  re'duction  m  ~  /'!•’ 


U,(t)e/,(X  -  tp)  -/( 


/'(  f(l  .)  Itl'i))  ^  s)  iioh) 


Te;  alio'v  fejr  a  re-eluctiejii  in  jiqe/dx)  —  /e  (.  we  rie-eel  to  she)w 


G6 


/'( -  -".'(Til  -}§')-  Hiis)  >  '*1  - 


A  small  value  of  say  0  <  <  grj.  implies  t  >  |-  A’ rider  the  structure  of  the 

algorithm,  this  rules  out  t  =  l/(9||p||^^),  a  value  of  t  necessary  for  the  previous 
convergence  and  complexity  jiroofs.  Therefore  in  this  case,  if  satisfies  it's  upper 
bound  of  convergence  and  polynomial  complexity  are  uncertain.  This  leads  to 
the  content  of  the  next  chapter,  where  analysis  yields  a  modification  to  the 
algorithm  which  allows  a  certain  quantifiable  reduction  in  u.'. 
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CHAPTER  VI:  MODIFYING  THE  PRIMAL  DUAL  ALGORITHM 

VI.  1  Modifying  the  Dual  V^ariable  Newton  Direction  and  Step  Size 

The  previous  chapter  showed  the  nature  of  the  Primal-Dual  Newton  direction 
in  the  dual  variables,  and  the  inherent  bias  and  uncertainty  in  its  effect  on 
reducing  juji7j(x)  — for  each  i  6  I.  This  chapter  develops  a  modification  to  the 
Newton  step  in  the  dual  variables  to  remove  this  bias  and  uncertainty.  We  begin 
with  the  simpler  case  of  an  arbitrary  linear  affine  5f,(x).  so  that 

+  tp)  -  p  = 

[gi(x)  +  Vgi(x)'^'tp)~fi. 

A  simple  way  to  modify  the  dual  New'ton  step  would  be  to  adjust  the  step  size  t 
(computed  in  the  manner  of  the  Primal  Variable  Algorithm)  using  information  at 
hand.  Let  this  adjustment  take  the  form  of  a  parameter  t-  multiplying  the  step 
size  t  in  u  (t).  Then  the  above  equation  becomes 
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Uv  1  - 


[^V5.(x)‘r,tp^  +  (^^  -  u.y^t  (5.(x)  +  V,9-(xrtp)- 


/I. 


=  (Ui.9,(x)  -  /0(1  -  r,t)  -  (u,^,{x)  -  /x)j^^^V5-(x)^tpy^t 


Now  lot  Tj  —  T-  =  - j- 


-.  With  this  value  we  obtain 


=>  u-(f;t)^.(x  +  tp)-// = 

(u,-5',(x)  -/^)(l  -  f^  t)  -(u-5-(x)  -/z)^^^V9-(x)'^tp^f;.t 
=  (Uj3j(x)  -/i)^l  -r-t^l  +-^Vg.(x)^tp^^ 

=  («i9,(x)  -/^)(1  -t). 


Thus,  multiplying  t  by  f)-  in  the  computation  of  the  new  value  of  u-  decreases 
|uj'5j(x)  —  /z  I  by  a  factor  of  (1  —  t)  at  each  inner  iteration.  (Recall  0  <  t  <  1.)  This 
has  a  powerful  implication.  Suppose  the  algorithm  is  at  a  point  near  the  trajectory. 
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that  is.  w^il  +  u  <  l|P||^v///i  <  and  n-g-{x)  - //  =  0. 

i  G  I.  Let  71  =  (1  -  9)g.  Then  u-g-(x)  -  /i  =  Ujfifj(x)  -  g -  dg/i  el.  If  each 
is  linear  affine  and  the  algorithm  uses  the  formula  Uj(f^  t),  after  any  given  iteration 
the  value  of  \i^g^[x)  —  g  is  the  same  for  all  i  G  I,  and  =  H  —  a;H.  where 
^’  =  1  —  Uj-5f^(x)//i,  i  G  I.  So  within  an  arbitrary  outer  iteration,  the  initial  common 
value  of  Uj5'j(x)  —  g  is  reduced  strictly  monotonicaily  at  each  inner  iteration. 
Furthermore,  if  any  inner  iteration  returns  a  step  size  t  =  1,  then  u-g-{x)  —  g  =  Q. 
i  G  I  at  each  subsequent  inner  iteration. 

This  result  applies  only  to  linear  affine  constraints,  so  the  next  logical  step 
would  be  to  attempt  a  .similar  result  for  the  more  general  concave  quadratic 
constraint.  We  begin  with  the  general  step  size  multiplier  r-  used  successfully 
before. 


5;(x) 


V^j(x)'^r-tp 


X  (3i(x)  + Vy-(x)'''tp+itp'^V“5,(x)tp)  -  g 


(6.1) 


r^‘.g»(x) 

9,(x) 


tp 


1 

5,(x) 


V5^,(x)’^’tp 


)) 


(6.3) 
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The  presence  of 


itp’ 


tp  in  (6.2)  not  multiplied  by  t-,  coupled  with  the 


possibility  that  V</j(x)^p  =  0  means  that  in  general  it  may  not  be  possible  to  zero 
out  (6.2)  and  (6.3)  with  the  correct  choice  of  r-.  But  the  components  containing 
in  (6.2)  and  (6.3)  can  be  made  to  sum  to  zero,  as  follows: 


Let  T- 


1  + 


.9/(x) 


and  the  following  equations  obtain. 


tp 


u,(''’j  t)gi(x  +  tp)  -/i  = 


/ 

(u-5-(x)  -  //)(  1  - 


1  V“5.(x) 


tp 


-Ui9,(x)  -2^P 


1  V  g-{x) 


.9i(x) 


tp 


=  -/0(1  -  t)  +  u,{itp'fV‘^5.(x)tp). 


The  use  of  t-  has  removed  .some  of  the  uncertainty  concerning  the  reduction  of 
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lu-iC'^x)  —  j,  but  a  negative  bias  remains:  u-Qtp*  V-5j(x)tp)  <  0.  However,  if  we 
define 


Uj(|tp''V-5,(x)tp) 


,  V“5.(x) 


tp 


5,(x 


we  obtain 


(u,(r- 1)  +  .i^)gi{x  +  tp)  -  p  =  u.(T.t)p^(x  +  tp)  -  p  +  J-g.(x  +  tp) 

=  -  p)(l  -  t)  +  U;(itp'V“5.(x)tp)-  U,-(itp'^'V2^.^X)tp) 


=  -p)(l  -t). 


The  use  of  r  ■  and  /3j  in  this  manner  removes  uncertainty  and  bias  concerning  the 
reduction  of  |uj-5’-(x)  —  p|  in  the  same  way  as  was  achieved  for  the  linear  affine 

n  — 

Pj(x);  and  when  V~g-{x)  =  0  for  all  x  (g-  is  linear  affine),  J,  =  0  and  this  technique 
reduces  to  the  method  developed  earlier  for  linear  affine  9,(x).  The  following 
definitions  will  allow  incorporation  of  this  modification  into  the  Primal-Dual 
Algorithm; 


u,(t,  T-,  ~  u-  -i- 


9ii^) 


9ii^) 


—  u. 


-t  -t-  J- 


or  in  vector-matrix  notation. 
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u(t,  T  ,  3  )=  u  +  ^  +  G  ^e/i  —  u  t  +  d 


where  3  is  an  m  x  1  vo?ctor.  the  entry  of  which  is  3- ,  and  r  is  an  m  x  m  matrix 
equal  to  diag(fj- ). 


VI.2  A  Modified  Primal-Dual  Algorithm 

The  Modified  Primal-Dual  Algorithm  (MPDA)  is  essentially  the  Primal 
Variable  Algorithm  with  the  addition  of  the  modified  Primal-Dual  Newton  step  in 
the  dual  variables.  Other  changes  have  been  made  so  that  the  requirement  of 
uj  <  1/62  is  met.  Consider  the  algorithm  at  the  end  of  an  outer  iteration  for  an 
arbitrary  ^  >  0.  meaning  the  tolerance  HPH^i^,; <  t  =  I'll  —  lj  has  been  met. 
Suppose  the  algorithm  has  also  achieved  Ujg^{x)  —  =  0,  i  €  L  This  implies  a;  =  0, 

M/p  =  H,  and  the  tolerance  thus  becomes  ||pll;^^ j ^  =  llPlIff  —  =  g-  point 

in  the  algorithm  that  g  is  decremented.  Let  Ji  denote  the  decremented  value,  so 
^  =  (1  —  for  some  0  G  (0,  1).  When  /x  is  so  decremented,  uj  no  longer  equals 
zero.  In  fact,  for  an  arbitrary  i  €  I, 


1  ( 

1  9 

-9 

n  1 

ll-^V 

i-e 

By  specifying  0  <  1/63,  we  ensure  that  when  reducing  p  to  /I  =  (l  —  0)/x  we  have 
u)  <  1/62  at  the  start  of  the  next  iteration.  Recall  that  with  the  modification  to 
the  dual  Newton  step,  uj  is  reduced  monotonically  during  the  course  of  an  outer 
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iteration,  and  once  a  step  size  t— 1  is  achieved  the  algorithm  maintains 


U'9j(x)  —  /i  =  0.  i  6  I  (and  therefore  uj  =  0)  for  the  remainder  of  that  outer  iteration 
until  the  tolerance  is  met  and  //  is  decremented.  If  the  situation  occurs  where  the 
tolerance  is  met  at  the  end  of  an  outer  iteration,  but  n-g-{x)  —  ^  0  for  some  i  G  I. 

the  algorithm  simply  takes  one  more  Primal-Dual  Newton  step  with  t  =  1, 
achieving  \i-g-{x)  —  /i  =  0,  i  G  I.  and  a;  =  0.  By  Corollary  IV. 3  this  maintains  the 
necessary  proximity  to  the  central  path,  g  is  then  decremented  to  /Z  =  (1  —  9)g 
resulting  in  u  =  1/62.  and  the  next  outer  iteration  begins. 


Input 

Reduction  factor  G  (0, 1/63] 

Tolerance  r  =  i'll  —  u,'.  .Assume  max 
^  i  G  I 


1 


alw’ays. 


T  is  an  accuracy  parameter. 

x*^  G  R*  is  given,  as  is  <  l/T,  and  u®  >  0,  i  G  I  such  that  |p(x®.u*^,//®)  <  t 

and  u?fif,(x'^)  - =  0,  i  G  I 


begin 


X  =  X®,  u  =  u®,  g  =  g^ 
while  g  >  T /2m  do 

begin  outer  iteration 


while  II  p  11^/^  <  r  do 

begin  inner  iteration 

A  =  arg  min  {$(x,n)  s.t.  (x  +  Ap)  G  IR*. 
0<A<l‘  ^ 

\\P\\M/g  ~  2wn|1  +  u;} 
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X  =  X  +  A  p 
u  =  u(  A  .  T  .  J  ) 
end  inner  iteration 


if  max 

1  e  I 


,  'i,f/-(x) 


>0 


begin  oik;  extra  step 
X  =  X  +  p 
u  =  n(  1.  T  .  J  ) 
(‘iid  one  extra  step 
/z  =  ( 1  -  6* )// 
end  outer  iteration 

end 


The  lemmas  and  convergence/complexity  theorems  apply  to  the  Modified  Primal- 
Dual  Algorithm  as  they  did  to  the  Primal  Variable  Algorithm.  The  bounds  on  the 
number  of  iterations  two  choices  of  0  are: 

^  =  =>  O^Vm|/n(e)|^  iterations 

0  =  0(1)  =»  0^mj/n(e)|^  iterations 

and  the  corresponding  computational  complexity  is: 
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0(n^>Jm|/n(e)|)  for  j 

O^n^m  |/n(e)|j  for  ^  =  0(1) 


As  stated  at  the  end  of  Chapter  3,  if  the  problem  P  has  only  linear  affine 
constraints,  then  setting  e  =  2  "  ^  allows  a  purification  algorithm  to  find  an  optimal 
vertex  with  no  increase  in  complexity  (see  Gonzaga(1992)). 

.\nother  significant  attribute  of  the  Modified  Primal-Dual  Algorithm  is 
convergence  of  the  Lagrangian  multipliers  to  an  optimal  vector  of  multipliers,  u*. 
This  results  from  the  way  in  which  the  MPDA  assumes  the  behavior  of  the 
standard  primal  variable  SUMT  at  the  end  of  each  outer  iteration.  Suppose  the 
algorithm  is  at  the  end  of  an  outer  iteration  for  an  arbitrary  /r  >  0.  Then 
n-g-{x)  =  /i,  i  G  I.  -Additional  modified  primal-dual  Newton  steps  will  coincide  with 
standard  SUMT  Newton  steps,  and  will  also  maintain  u,-^j(x)  —  g,  i  el.  Thus,  the 
MPDA  inherits  the  Lagrangian  multiplier  convergence  property  of  SUMT,  as 
proved  in  Theorems  25  and  26  of  Fiacco  and  McCormick(1968,  1990). 
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CHAPTER  VII:  LINEAR  PROGRAMMING  WITH  THE  MODIFIED 
PRIMAL-DUAL  ALGORITHM 

\TI.1  A  Review  of  Aiistreieher's  Results 

In  previous  chapters  we  have  developed  a  modification  to  the  Primal-Dual 
Newton  step  which  allowed  for  a  guaranteed  reduction  in  the  deviation  from 
perturbed  complementary  slackness  with  each  Newton  step.  That  reduction  in 
turn  allowed  the  application  of  the  den  Hertog,  Roos,  and  Terlaky(1990)  analytical 
techniques  to  prove  convergence  and  polynomial  complexity  for  the  Modified 
Primal-Dual  Algorithm.  This  chapter  will  show  that  the  modificationt  also  allows 
us  to  adapt  Anstreicher’s  results  for  quadratic  programming  with  the  logarithmic 
barrier  function  to  linear  programming  with  the  Modified  Primal-Dual  Algorithm. 
Thus,  we  have  another  way  to  show  polynomial  complexity  for  the  Modified 
Primal-Dual  Algorithm,  albeit  for  a  less  general  class  of  problems.  (An  attempt 
was  made  to  adapt  Anstreicher's  results  directly  to  the  Primal  Variable  Algorithm 
for  small  w,  as  was  done  with  the  work  of  den  Hertog,  Roos,  and  Terlaky.  Because 
of  the  relatively  tight  nature  of  the  inequalities  in  several  of  the  lemmas,  however. 
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this  has  not  been  successful.) 

We  begin  by  reviewing  Anstreicher's  convergence  and  complexity  proofs  for 
the  traditional  SUMT  logarithmic  barrier  algorithm  applied  to  a  quadratic 
program.  Such  a  review  allows  the  reader  to  see  a  significantly  different  approach 
to  proving  polynoniality,  one  that  makes  important  use  of  the  linear  algebra 
associated  with  the  linear  affine  constraints  of  the  quadratic  program.  W'e  then 
show  that  the  Modified  Primal-Dual  Algorithm,  when  applied  to  a  linear  program, 
can  be  analyzed  for  convergence  and  complexity  in  the  same  manner  as 
Anstreicher's  treatment  of  SUMT.  The  lemmas  and  theorems  from 
.A.nstreicher(T990)  are  presented  in  this  chapter  without  proof.  For  the  interested 
reader,  the  proofs  are  in  .Appendix  1.  precisely  as  developed  in  Anstreicher(1990) 
and  .Anstreicher  et  al(1990)  with  only  notational  changes  to  conform  to  this 
dissertation. 

Let  us  define  the  standard  quadratic  program  (QP)  in  using  notation 
similar  to  Anstreicher: 

QP;  min  f(x)  =  6^x  +  |x^Qx 
s.t.  .Ax  —  c  >  0 

for  b  €  E”,  c  G  E”,  .A  an  m  x  n  matrix,  and  Q  an  n  x  n  symmetric  positive  semi- 
definite  matrix.  We  assume  the  data  is  integer  with  bit  length  L.  and  as  before  the 
feasible  region  R  is  bounded  with  a  non-empty  interior  R*.  This  implies  that 
m  >  n  and  the  columns  of  A  are  linearly  independent.  The  WMlfe  dual  of  QP  is 

QD:  max  b^x -b  ^x^Qx  —  u^(Ax  —  c) 
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s.t.  b  +  Qx  —  A  ‘n  =  0 


i  u  >0 

I 

The  dual  can  bo  simplified,  b  +  Qx  —  A  ‘  u  =  0  =>  u  ‘  A  =  b  ^  +  x 


b  ^x  +  |x  ^  Qx  —  u  ^Ax  -I- 11  ’ c  =  b  ‘x  +  ix ^Qx  -  (b^’  +  x ^  Q)x  +  c 
=  —  |x^  Qx  +  c 


Thus  we  can  write  the  dual  problem  as 

QD:  max  d(x,u)  =  —  ^xHjx  +  c^u 
s.t.  b  +  Ox  —  A  =  0 
u  >  0 


Lemma  VII.T.  Let  x  be  feasible  in  QP  and  (x,  u)  be  feasible  in  QD.  Then 
f(x)  —  d(x,  u)  >  0, 

Let  denote  the  i^^  row  of  A.  We  define  the  logarithmic  barrier  function 
$(x,  for  x  €  (R*  and  //  >  0  as  before. 

^  -  E  in{a-'^x  -  C-) 

1=1 

The  gradient  and  Hessian  of  are 
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V$(x./i) 


Qx  +  6  >n  a- 
^  ,=l(«rx-c-) 


The  li  X  n  matrix  H  will  again  denote  the  Hessian  of  Phi.  By  Lemma  III. 2.  H  is 
positive  definite.  ‘I>(.x./i)  is  continuous  on  IR*  and  goes  to  +  oo  at  the  boundary  of 
R.  so  the  compactness  of  R  guarantees  $(x,/i)  has  a  minimizer  x(/x)  6  IR*. 

For  an  arbitrary  /i  >  0.  Anstreicher  employs  three  different  measures  of  the 
distance  of  a  point  x  6  R*  from  the  central  path  x(p).  Some  notational  definitions 
are  necessary:  s  =  Ax  —  c.  S  =  diag(s),  s  "■  *  =  (l/s[,  l/so,  .  .  ..  1/s^)^  =  S  “  ^e. 

The  primal  variable  Xewton  step  is  pj.  =  —  H  ~  ^  V$(x.p).  and  pg  =  Apj.  is  the  step 
in  s  resulting  from  the  Newton  step  p^  in  x.  The  three  measures  are; 


(1)  d(x.p) 


mm 


A^u  =  6  +  Qx 


^Su 


(2)  i|S-‘pg 


(3)  II Pi  II// 

To  better  understand  these  measures,  it  may  be  helpful  to  apply  them  to  the 
familiar  problem  P  from  previous  chapters:  ]|pSu  -  e||  is  the  same  as  \\n^g  —  p||,  (3) 
is  the  norm  of  the  Newton  step  in  the  primal  variables  with  respect  to  the  Hessian 
of  $(x,//),  and  (2)  is  the  same  as  (3)  but  ignoring  the  contribution  of  the  objective 
function  f(x). 
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Lemma  VII.2:  Let  x  G  R*  and  /x  >  0  be  arbitrary.  Then 

||S“‘ps|P<||Px||//  <<5(x./i)“. 

Lemma  VII.S:  Given  x  G  R*.  suppose  ||S  ~  <  1.  Let  x  =  x  +  p_p.  s'  =  s  +  p^, 

and  S  =  diag(^  ).  Let  H  =  V“$(x  ,//),  pj^  =  —  H  ~  ^  V$(x  .//),  and  j>j  =  Apj  . 
Then  x  G  R*  and 

II  S  ~  ‘ps  ||<llPx  !|'^  <  <^(x  .p)  <||S“  Vsp  <||Pj;|I7/  <  6ix.fl)- 

ii  ir 

Lemma  VTT  4:  If  x  G  R"^  and  i|Px||//  ^  T  then  ^(x.fi)  —  4>(x(/z)./i/)  <  — — . 

l-||Px||ff 

Lemma  VIL5:  If  x  G  R*  and  ||Px||^  <  !•  then 

|f(x)-f(x(rt)i<iiy^4’i^jM/LU 

The  algorithm  used  by  Anstreicher  and  Anstreicher  et  al  is  similar  to  that  of 
den  Hertog,  Roos,  and  Terlaky(1990).  It  begins  with  a  feasible  interior  point  x® 
and  an  initial  barrier  parameter  value  /x®  =  From  point  x^,  the  algorithm 

conducts  linesearches  along  a  sequence  of  Newton  steps  until  it  reaches  a  point 
where  ||  Px  ||//  <  t.  This  point  becomes  x^  "*■  * ,  is  decremented  to 

—  9)fi^  for  some  fixed  B  G(0, 1),  and  the  process  starts  over  again.  The 
algorithm  terminates  when  ||Px||^  S  f  =  2  from  which  the 
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optimal  solution  may  be  found  in  O(m^)  operations.  Convergence  and  complexity 
of  the  algorithm  arc  proven  by  Aiistreicher  with  the  following  two  theorems. 


Theorem  VTI.l: 


Let  z*  denote  the  minimum  value  of  f{x)  for  problem  QP. 


After 


0(L/^)  outer  iterations,  the  algorithm  obtains  ^  such  that  m/i^’  <  2  and 

a  point  X  =  x^'  ^  such  that  f(xi  —  z*  <  ^  ^ 


The  proof  of  Theorem  VII. 2  requires  the  following  lemma  proving  a  guaranteed 
reduction  in  $(  •  ,/i]  along  the  Xewton  direction.  It  is  similar  to  the  den  Hertog, 
Roos.  and  Terlaky  Taylor  series  expansion  result  used  in  the  proof  of  Lemma  IV.  1. 
In  that  case,  the  proof  was  quite  long  and  detailed,  and  was  omitted.  .Anstreicher’s 
proof  for  this  lemma  is  somewhat  shorter,  and  is  included  in  Appendix  1. 

Lemma  VTI.6;  Let  A  =  -4 — — .  Then 

- 1+llPxll// 

$(x,/i)  -$(x  +  A  p^,/x)  >jjp^||^-/n{l  +||p^||^) 

Corollary  VII. 6:  Let  r  <  ||Pa:||//  <  Then  a  pure  Newton  step  (A  =  1)  obtains 
$(x,/i)  -  $(x  +  PxiM)  >  ^  —  r)  +  T  >  0 


We  can  now  proceed  to  Theorem  VII.2,  which  bounds  the  number  of  inner 


iterations  per  outer  iteration. 


Theorem  VII. 2:  Let  A  =  t  —  /n(  1  +  r),  E/-  =  and  N  denote 

the  number  of  inner  iterations  required  for  an  arbitrary  outer  iteratation.  Then  N 
is  bounded  as: 


N 


< 


{6m  +  2E ^^^m^  +  Ejj,  j 


Combining  Theorem  VH.l  and  VH.2  gives  the  complexity  result  we  seek:  the  total 
number  of  Newton  iterations  required  by  the  algorithm  to  solve  problem  QP  is 
bounded  by 


^  ^(9m  i- 2EyVm  ^  e<{,^ 

^  TZTe  +T'/ 


0(L). 


The  choice  of  6  has  much  the  same  influence  as  in  the  results  of  den  Hertog,  Roos, 
and  Terlaky: 


6  =  =$•  0(^fm  L)  iterations 


6  —  0(1)  =>  O(mL)  iterations. 


In  Corollary  VII. 6  it  was  shown  that  a  pure  Newton  step  would  achieve  a 
guaranteed  reduction  in  $(  ■  ,/i)  if  r  <  |[Pxl'iy  <  To  ensure  that  ||Px||^  <  ^ 
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throughout  the  algorithm,  note  that  Theorem  VII. 2  implies  that 


is  the  maximum  reduction  possible  for  any  step  size,  while  Lemma  VII. 6  shows  a 
reduction  of  0.5  —  /n(1.5)  >  .0945  is  guaranteed  with  the  correct  choice  of  step  size 
for  II  p^  11^  >  Thus  it  suffices  to  restrict  the  maximum  possible  reduction  to  less 
than  .0945  to  ensure  that  ||Pj||//  <  This  may  be  done  by  choosing  6  such  that 

+  2EyN[mj  +  E<j,  <  .0945.  (7.1) 

For  small  values  of  r  both  Ej  and  E^  approach  zero,  so  9  such  that  9‘^m  <  .0945,  or 
9  <  .3/^Jm  will  suffice.  This  result  will  be  important  in  later  analysis. 

This  concludes  the  review  of  Anstreicher's  complexity  results  for  QP;  the 
remaining  portion  of  this  chapter  will  show  that  these  results  may  be  applied  to 
linear  programming  with  the  Modified  Primal- Dual  Algorithm  to  achieve  the  same 
complexity  bound. 

VII. 2  Another  Polynomial  Complexity  Result  for  the  Modified  Primal  Dual 
Algorithm 

To  remain  consistant  with  the  notation  used  in  the  development  of  the 
Modified  Primal-Dual  Algorithm,  let  us  redefine  the  linear  program  (LP)  in  a 
manner  identical  to  the  definition  of  problem  P  in  Chapter  III,  with  the  additional 
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restriction  that  the  constraints  be  linear  affine: 

LP:  min  b'^x 

s.t.  gj(x)  >  0  and  linear  affine,  i  E  I 
X  6  E” 

All  restrictions  on  the  feasible  region  for  P  apply  to  LP.  This  LP  is  the  problem 
we  will  address  throughout  the  rest  of  this  chapter. 


In  Chapter  III  this  thesis  introduced  the  matrix  M  and  its  relationship  to  the 
Hessian  of  the  logarithmic  barrier  function.  Of  particular  interest  is 


where  Uj  >  0,  i  G  1.  As  before,  let  H  denote  the  Hessian  of  $(  •  ,/r).  It  is  clear  that 
if  Uj-5^(x)  —  /r  =  0  for  all  i  6  I,  then  H(x,/i)  =  ^M(x.^). 

Lemma  VIL7:  Let  u^gj(x)  —  fi  =  0  for  all  i  €  I,  and  let  ^  G  (0, 1)  be  given.  If 
J[  =  (l  —  then 


iM(x,/x  )  H(x,^) 


Proof:  From  (3.7)  we  have 


rVg,(x)Vff.-(x)'" 
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Using  Uj5j(x)  =  i^i  and  /z  =  ( 1  —  9)^,  we  get 


-  rn 

^M(x.^)=  £■ 


/i 


Vg,(x)Vff,(x)'^ 

5/(x)‘^ 


=  (rr?)HM 


which  completes  the  proof.  □ 


A  Linear  Programming  Version  of  the  MPDA.  The  algorithm  begins  with  a 
feasible  interior  point  x®.  an  initial  barrier  parameter  value  p®  =  >  0  such 

that  u^g-{x*^)  —  =  0  for  all  i  E  I.  and  ||P||,\//^  —  ||P||/y  ^  Let  the  outer 

iteration  be  arbitrary.  Starting  from  point  x^,  the  algorithm  takes  a  series  of  pure 
Primal-Dual  Newton  steps  in  the  primal  variables,  that  is, 


At  each  such  step,  each  Uj,  i  €  I  is  updated  to  u,(l,f^-  ,0).  These  inner  iterations 

continue  until  the  algorithm  reaches  a  point  where  HPill^/  <  t,  say  the  iteration. 

Then  x^  ^  =  x-,  is  decremented  to  ^  =  (l  —  9)fi^  for  some  fixed  ^  6  (0, 1) 
j 

satisfying  (7.1),  and  the  process  starts  over  again.  The  algorithm  terminates  when 
II Px II//  ^  ^  and  xnn^  <€  =  2~^^^\ 

The  initial  Primal-Dual  Newton  step  in  the  primal  variables  will  be  identical  in 
direction  to  that  generated  by  the  logarithmic  barrier  function,  but  will  be  shorter 
by  a  factor  of  (1  —  6).  Therefore,  the  reduction  A  in  $(  •  ,/i)  from  this  first  iteration 
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will  be 


A  >  (1  -^)r2  +  /n(l  -(1  -^)r)+(l  -0)r.  (7.2) 

Since  T  is  small  and  positive,  the  right  hand  side  of  (7.2)  is  strictly  positive  for 
0  =  0,  equals  zero  for  0  =  1.  and  is  strictly  monotonically  decreasing  with  respect  to 
0  on  (0, 1).  Thus,  A  >  0  for  0  6  (0,1).  After  the  initial  pure  Primal-Dual  Newton 
step  in  the  primal  variables.  M/p  =  H  and  the  exact  results  of  Anstreicher  are 
obtained.  (It  is  also  permissible  to  conduct  linesearches  along  the  Primal-Dual 
Newton  directions  after  the  first  inner  iteration's  pure  Primal-Dual  Newton  step. 
The  decrease  in  <5(  • .//)  will  then  be  guaranteed  by  Lemma  VII.6  instead  of 
Corollary  "VII.G.)  Therefore,  the  complexity  results  proved  by  Theorems  VII.l  and 
VII. 2  for  the  logarithmic  barrier  function  applied  to  QP  apply  for  this  version  of 
the  Modified  Primal-Dual  Algorithm  applied  to  LP. 


87 


CHAPTER  VIII:  SATISFYING  THE  INITIAL  CONDITION  REQUIREMENTS 

We  begin  this  c  hapter  by  stating  the  general  problem  to  be  addressed  and  key 
notation  used. 

min  f(x) 
xeE’^ 

s.t.  5,(x)  >  0,  i  €  I  =  (1,  2,  .  .  m) 

where  —  f(x)  and  5rj(x),  i  6  (1,  2,  .  .  m)  are  concave  quadratic  functions. 

Standard  Notation  and  Definitions 

R  =  {x  G  E"  I  ^j(x)  >  0,  i  G  1} 

R*  =  (x  G  E"  I  5,(x)  >  0,  i  G  1} 

i=l 
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Vf(x)  ^Vg,(x) 


Hfx.;,)  ^  V2$(x,,,)  =  g^9,(x)^9.(x)’'  V^s.lx) 


itl  '7,(x) 


2 


9ii^} 


(Note:  when  f(x)  is  linear.  V"f(x)  vanishes.) 


lit 

M(x.u./i)  = 


V,9,(x)V.9,(x)'"  V2g,(x 


i  =  i 


g-(xf  9rM 


p(x,u./i)  =  -M' 


b'‘  -  /i2  V5,(x)/^,(x) 
>=1 


=  -(1m)-'V4.(x,^) 


When  clarity  is  not  affected,  let  p  =  p(x,u,^),  H  =  H(x./i),  and  M  =  M{x,u,/z). 


IIPII^  =(p'^Hp)^/2 


For  some  zeE^  let  the  region  Rj  =  {xcE”  |  f(x)  <  z,  x  e  R}  be  bounded  with  a  non¬ 
empty  interior.  Then  for  some  positive  integer  r  we  define  the  analyt^ic  center  of 
R^  as  the  point  in  R^  whi..h  solves 

max  (z -f(x))'’n  5j(x). 

i=l 
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This  dissertation  has  presented  two  implementations  of  the  Modified  Primal- 
Dual  Algorithm  (MPDA),  demonstrating  two  ways  to  show  polynomial  complexity 
for  the  algorithm.  In  each  of  these  implementations,  the  analysis  assumes  a  fairly 
restrictive  initial  condition;  nameh  .  a  known  initial  point  (x^.u®)  and  barrier 
parameter  p®  >  0  such  that  the  following  conditions  hold: 


x"  G 


(8.1) 


(8.2) 


(8.3! 


Recall  that  if  u!-g^{x^)  —  /r  =  0,  i  G  I,  then  (8.2)  is  equiv-.^lent  to  ||p(x*^.u*^,/r®)  <  r. 

.Also.  (8.3)  applies  only  for  problem  LP,  that  is.  when  the  objective  function  is 
linear  and  the  constraints  linear  altine.  In  practice,  one  may  not  receive  such  a 
conveniently  structured  prof  ie:.  Fortunately,  it  is  a  straightforwara  task  to  meet 
these  conditions  given  an  arbitrary  starting  point  x  G  E’^.  Our  first  concern  is 
finding  a  feasible  value  of  x. 

VIII. 1  Feasibility 


The  following  method  for  achieving  feasibility  is  a  version  of  the  Two  Phase 


method  often  used  in  linear  programming.  For  a  detailed  discussion  of  that 


application,  see  Bazaraa.  Jarvis,  and  SlieraJi(  1990).  Given  x,  let  V  denote  the  set 
of  vndicies  of  constraints  that  are  non-positive  for  x.  Xow  define  a  set  of  artificial 
variables  Oj,  j  6  and  a  variable  z  £  E*.  Let  =  -  (7j(x)  +  1.  j  £ 
z^  =  ma.x{aj)  +  1.  and  foriniilate  the  following  problem  FP. 

FP:  min  z 

s.t.  .9j(x)  >  0,  i  £  I  -  T 

+  >  0.  j  £  r 

z  >  Gj.  J  £  r 
z  >  - 1 

By  construction.  FP  has  a  solution  and  (x.  a®,  z®)  is  strictly  feasible  for  FP. 
Furthermore.  FP  has  a  linear  objective  function  and  concave  quadratic  constraints, 
so  the  algorithm  of  den  Hertog.  Roos,  and  Terlaky  can  solve  FP  in  polynomial 
time.  The  number  of  constraints  remains  0(m),  and  if  we  assume  m  =  0(n;,  the 
computational  complexity  is  the  same  as  for  problem  P.  In  reality,  the  algorithm 
may  be  terminated  as  soon  it  obtains  a  negative  value  of  z.  The  associated  value  of 
x  is  strictly  feasible  for  P  and  may  be  denoted  x  .  If  the  minimum  z  is  non¬ 
negative.  then  P  is  infeasible.  With  Sc  £  IR*,  we  have  satisfied  (8.1),  so  we  now 
proceed  to  satisfying  (8.2)  and  (8.3),  known  as  centering. 

Another  approach  to  finding  an  interior  point  is  found  in  Fiacco  and 
McCormick(1968,  1990).  Let  x  £  E”  be  arbitrary,  and  define  the  sets 
S  =  {i  I  5'j(x)  <  0,  i  £  1}  and  T  =  {i  1  5j(x)  >  0,  i  £  I}.  Now  minimize 

U(x,//^')  =  -  X]  + 

ieS  i^T 


91 


for  /i^  jO  as  k  -►  -f  oo.  When  in  the  course  of  this  minimization  a  new  constraint  is 
strictly  satisfied,  its  index  is  moved  from  S  to  T,  is  redefined  accordingly, 

and  the  process  is  begun  again.  This  continues  until  all  constraints  are  strictly 
satisfied  or  it  is  shown  no  interior  point  exists. 

VHI.2  Centering 

A  straightforward  method  for  initial  centering  is  found  in  McCormick(1991b), 
based  on  the  work  of  Huard(19G7)  and  Fiacco(1979).  Given  that  "x  G  R*  (possibly 
from  the  application  of  the  feasibility  methodology  above),  define  V®  =  {x  G  E”  | 
f(x)  <  f(  ){  )}  n  IR*.  One  then  finds  an  approximate  .solution  to  the  problem 

max  [ -f(x) +  f(  X  )]n5i  (x).  (8.4) 

X  G  «=l 

This  is  equivalent  to  finding  an  approximate  solution  to 

min  -/n[-f(x) +f(  X  )]- 2/n(g  (x)).  (8.5) 

xGV"  »=1  ^  ’ 

Let  x^  approximately  minimize  (8.5).  Then  we  have 

Vf(x^)  m  Vg,(x*^) 

[  -  f(x«)  +  f(  X  )]  .^1 

approximately  equal  to  zero.  Now  let  =  f(  Sc  )  —  f(x®).  and  we  have  (x®,//®)  that 
approximately  solves 
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min 
X  e  R* 

With  the  (familiar)  definition  of  u®  =  and  given  a  good  enough 

approximate  solution,  we  have  (x^,  u*^,  fx^)  which  will  satisfy  (8.2). 

Although  this  two  step  approach  is  a  standard  one  and  very  workable,  we  now  turn 
to  a  method  which  combines  the  feasibility  and  centering  operations  and  makes  use 
of  a  polynomial-time  analytical  centering  algorithm.  The  method  is  based  on  Jarre 
(1991),  which  demonstrated  the  polynomality  of  a  method  of  analytical  centers 
when  applied  to  our  problem  P.  The  presentation  here  will  noc  include  detail  on 
the  polynomality  of  the  centering  algorithm,  which  is  beyond  the  scope  of  this 
research- -readers  interested  in  such  detail  are  better  served  by  reading  the  original 
paper.  Rather,  we  show  a  scheme  adapted  from  Jarre  (1991)  which  allows  the 
application  of  his  algorithm  to  finding  (x®,  u®,  satisfying  the  initial  conditions 
of  the  MPDA. 


Combined  Feasibility  and  Centerine 


We  liiake  the  following  definitions: 


d-  =  g-{0)  -  1,  i  G  I 

=  9ii^)  ~  id-,  i  G  I 
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rn 

^m+l^  -  E  ^9,(0) 

2  =  1 


This  gives  us 


9,(0,1)  =  Life  (1,2 . m+1) 

h  ~  ■ 


(8.6) 


Now  let 


1R(7,i/)  =  (x  €  E”  I  9i(x,7)  >  0,  i  €  I,  9m  +  ^  0}. 

for  some  fixed  pair  (7,1^)  such  that  7  6  E^,  u  <00.  Since  R  is  a  bounded  convex 
polytope  and  1R(7,!^)  has  as  constraints  perturbations  of  all  the  constraints  forming 
(R,  ^(7,1^)  must  be  bounded  (Fiacco  and  McCormick(1956,  1990),  Theorem  24). 

The  point  x  =  0  is  an  element  of  R(l,l)  so  R(l,l)  is  non-empty.  In  fact,  (8.6) 
shows  that  x=0  is  the  analytical  center  of  R(l,l).  We  now  make  two  successive 
applications  of  the  analytical  centering  algorithm,  the  first  with  a  strictly 
monotonically  increasing  sequence  {i/^}too  with  1/^  =  1,  and  the  second  with  a 
strictly  monotonically  decreasing  sequence  with  7*^  =  1.  The  first 

application  of  the  algorithm  “backs  off”  the  m+1^^  constraint  until  it  no  longer  has 
any  significant  effect  on  the  analytic  center,  and  the  second  application  decreases 
the  perturbations  to  the  first  m  constraints  until  R  closely  approximates  R,  and 


9' 


thus  the  resulting  close  approximation  of  the  analytical  center  of  [R  is  a  close 
aproximation  of  the  analytical  center  of  K.  In  this  way.  two  applications  of  a 
polynomial  analytical  centering  algorithm  yield  a  a  close  aproximation  of  the 
analytical  center  of  IR,  denoted  x®,  which  implies 


is  nearly  zero.  All  that  remains  is  to  then  select  a  large  value  of  p®  such  that 
^0  ^  so  that  for  problem  P. 


remains  close  to  zero,  and  the  usual  definition  of  u®  =  yields  (x*^,  u*^,  ffi) 

satisfying  the  initial  conditions  in  (8.1)  through  (8.3). 
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CHAPTER  IX:  SUMMARY  AND  RECOMMENDATIONS  FOR  FUTURE 
RESEARCH 

IX. 1  Summary  of  Results 

We  begin  this  summary  by  restating  the  problem  to  be  solved  along  with 
important  notational  conventions. 

Standard  Notation  and  Definitions 

General  Constrained  Optimization  Problem: 

min  f(x) 
leE^ 

s.t.  9-(x)  >  0,  i  G  I  =  (1,  2,  .  .  m) 

R  =  {x  G  E”  I  ^j(x)  >  0,  i  G  1} 


96 


R*  -  {x  e  E"  I  g^ix)  >  0,  i  6  1} 


$(x,//)  =  for  X  6  0?*  and  ^  >  0. 

«=1 


H  =  V^(x.g)  = 


V"'f(x) 


,  ^V5-(x)V5-(x)'^ 
1=1  Sii^) 


V^gi(x) 

ffi(^) 


(Note:  when  f(x)  is  linear,  V^f(x)  vanishes.) 


L(x,u)  =  f(x)  -  2  Uj^j(x) 
i=l 


VL(x,u)  =  Vf(x)  -  2  u,V^-(x)  =  0,  u,  >  0,  i  e  I 
i=l 


m 

t=i 


V^.{x)V(7.(x)'r 

V^S.lx)' 

ffi(x) 

for  X  G  R*  and  Uj  >  0,  i  G  I. 


P 


-M 


1  -  /i  2  Vg,(x)/gi(x) 

f=l 


-(iM)-*V$(x,/i) 
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u 


max 
i  G  I 


^  9 


-^IIPII//  ^  hi  1(1  ^  'll  IIP II//  for  u;  <  1  and  f(x)  linear 


This  dissertation  began  with  a  review  of  concepts  fundamental  to  constradned 
optimization  with  interior  point  methods  and  current  research  in  the  field.  That 
was  followed  by  a  detailed  description  of  the  Primal-Dual  Algorithm  developed  by 
McCormick(1991b).  Recall  the  basic  idea  of  that  algorithm  is  to  use  Newton’s 
method  with  a  step  size  line  search  to  find  solutions  to  the  following  set  of 
equations: 


Va.L(x,u)  =  0 


-  /X*  =  0,  i  G  I 


where  x  G  R*  and  u,  >  0,  i  G  I.  These  equations  are  solved  for  a  sequence  of 
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positive  fj.^  such  that  as  k->-f  oo.  The  first  major  research  result,  in 

Chapter  IV,  was  showing  the  similarity  of  p,  the  primal  variable  Newton  direction 
from  the  Primal-Dual  Algorithm,  with  the  Newton  direction  —  H  “  ^  V$(x/i)  for  the 
standard  SUMT  with  the  logarithmic  barrier  function.  This  similarity  was 
exploited  by  constructing  a  Primal  Variable  Algorithm  using  the  search  direction  p, 
and  limiting  the  divergence  of  p  from  —  H  ~  ^  V$(x//)  through  the  measure  u.  It 
was  shown  that  by  keeping  w  small  {u  <  1/62),  the  SUMT  polynomiality  results  of 
den  Hertog,  Roos.  and  Terlaky(1990)  could  be  used  as  a  template  to  construct  a 
convergence  and  polynomial  complexity  proof  for  the  Primal  Variable  Algorithm. 

In  the  Primal  Variable  Algorithm  analysis,  ui  <  1/62  was  specified.  Ensuring 
that  uj  meets  that  criterion  is  another  matter.  The  degree  to  which 
u,5f,(x)  —  =  0,  i  g  I  is  violated  during  the  course  of  the  PrimaJ-Dual  Algorithm  is 

very  dependent  on  the  way  the  Uj  Vciriables  change.  Thus,  Chapter  V  contained 
research  into  the  performance  of  the  ducd  variable  Newtcn  direction  for  an 
arbitrary  u,,  particularly  when  using  a  step  size  generated  by  the  Primal  Variable 
Algorithm.  The  result  was  two  cases,  determined  by  the  sign  of  Uj^,(x)  —  Both 
cases  revealed  that  the  maintenance  of  u;  <  1/62  is  not  certain  when  movement  in 
the  dual  variables  Uj,  i  G  I,  according  to  the  dual  Newton  direction,  is  appended  to 
the  Primal- Variable  Algorithm.  Thus,  further  study  was  required  to  establish 
polynomiality  for  the  Primal-Dual  Algorithm. 

In  Chapter  VI,  research  centered  on  developing  a  modification  to  the  dual 
variable  Newton  step  for  an  arbitrary  u,  to  guarantee  u;  <  1/62.  First,  the  simpler 
case  of  a  linear  affine  g',(x)  lead  to  the  discovery  that  a  simple  step  size  change 
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using  existing  information  could  guarantee  a  monotone  decrease  in  |uj5j(x)  — 
with  each  Newton  step.  Extending  this  result  to  the  more  general  quadratic  5j(x) 
required  both  a  step  size  change  and  an  additive  term  to  the  dual  Newton  step, 
again  with  existing  information.  With  a  monotone  decrease  in  juj5j(x)  —  at 
each  Newton  step  now  ensured  for  i  €  I,  movement  in  the  dual  variables  according 
to  the  modified  to  the  dual  Newton  step  was  added  to  the  Primal  V^ariable 
Algorithm,  resulting  in  the  Modified  Primal-Dual  Algorithm  (MPDA).  The 
MPDA  is  essentially  the  original  Primal-Dual  Algorithm  with  the  modified  dual 
Newton  step,  using  $(x,/i)  as  a  merit  function.  The  modified  dual  Newton  step 
and  small  decrements  to  ^  guarantee  ui  <  1/62,  so  convergence  and  polynomiality 
follow  from  the  Primal  V'ariable  Algorithm. 

The  effect  of  the  modification  to  the  dual  Newton  step  had  further 
consequences.  In  Chapter  VII,  analysis  revealed  that  when  the  MPDA  is  applied 
to  a  linear  program  (LP),  the  polynomiality  results  of  Anstreicher(1990)  apply 
directly.  Anstreicher's  work  is  remarkable  for  showing  polynomiality  of  the  original 
SUMT  of  Fiacco  and  McCormick  for  quadratic  programming  (QP)  when  initial 
conditions  and  parameter  values  are  correctly  specified.  When  a  linear  objective 
function  is  specified  to  conform  to  our  MPDA  analysis,  QP  becomes  LP,  and  it  is 
shown  that  Anstreicher’s  proofs  extend  to  the  MPDA.  Hence,  we  have  another 
convergence  and  polynomial  complexity  result,  although  for  a  less  general  class  of 
problem. 

The  research  of  this  dissertation  concludes  in  Chapter  VIII,  which  addressed 
the  problem  of  finding  an  initial  point  (x®,u®,/x®)  which  meets  the  fairly  restrictive 
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initial  conditions  used  by  both  den  Hertog,  Roos,  and  Terlaky(1990)  and 
Anstreicher(1990).  The  conditions  are  feasibility,  namely  G  [R*.  >  0,  and 

/>0;  and  centering,  meaning  ||  p(x  is  small.  Techniques  to  achieve 

feasibility  are  covered  by  many  authors-feasibility  results  from  Bazaara.  Jarvis, 
and  Sherali(1990)  and  Fiacco(1979)  apply,  and  were  covered  in  Chapter  VIII.  The 
chapter  also  presented  a  centering  method  from  McCormick(  1991b),  based  on  the 
work  of  Huard(1967)  and  Fiacco(1979),  and  concluded  with  a  combined  feasibility 
and  centering  technique  adapted  from  Jarre(1991). 


IX. 2  Future  Research 

It  is  possible  that  the  results  developed  here  for  convex  programming  with 
quadratic  constraints  could  be  extended  to  smooth  convex  programming.  We 
define  smooth  convex  programming  by  way  of  a  relative  Lipschitz  condition  set 
forth  in  Jarre(1991).  Let  x  and  x  +  h  in  R*  and  z  G  E”  be  arbitrary.  Then  we  say 
the  relative  Lipschitz  condition  is  satisfied  if  for  some  M  >  0,  the  Hessian  of  each 
5'j(x)  satisfies  the  following: 

+  h)  -  V2y.(x)]  z|  <  M||h||^(^jz’^V25,(x)z 


where 


H(x) 


m  Vg.|;x)Vg-(x)^ 

h 
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The  convergence  and  complexity  results  for  the  logarithmic  barrier  function  in  den 
Hertog,  Roos,  and  Terlaky(1990)  extend  to  smooth  convex  programs,  with  M 
entering  into  the  complexity  polynomial.  It  should  be  straightforward  to  make  the 
same  extension  to  the  Modified  Primal-Dual  Algorithm. 

Another  research  area  is  the  choice  of  merit  function.  For  the  general 
constrained  optimization  problem 

min  f(x) 

1  e  E^' 

s.t.  g'{x)  >  0,  i  e  I. 

there  are  two  candidate  functions  which  incorporate  primal  and  dual  variables 
mentioned  in  McCormick(1991a), 

T7i  m  I-  -i 

MFl(x,u,/i)  =  f(x)  -  mE  +  EKifjlx)  -  /i/u(u  -5f  (x))l,  and 

i=l  i=V-  *  J 

MF2(x,u,/z)  =  VL(x,u)'^M  ~  ^ VL(x,u)  +  £  ^u-^-  -  gf, 

i=l 


for  x  €  R*,  u  >  0,  and  /x  >  0,  where  the  matrix  M  is  defined  as  in  earlier  chapters: 


m 


i=l 


Vgi(x)Vg,(x)'^  v2g.(x) 


9i{^y 


9M) 


for  X  G 


and  Uj  >  0.  i  €  I. 
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MFl(x,u,/i)  has  the  interesting  property  in  that  a  stationary  point  with  respect  to  x 
and  u  for  MFl(x,u,/i)  satisfies  the  equations  motivating  the  Primal-Dual  Algorithm 
(McCormick(  1991b)); 

VL(x,u)  =  0 

Uf5e(x)  -  /i  =  0.  i  e  I. 

Successful  use  of  such  a  merit  function  might  reduce  or  eliminate  the  need  for  a 
modification  to  the  Primal-Dual  Algorithm  to  show  polynomiality. 

Finally,  only  limited  mention  has  been  made  in  this  dissertation  to 
convergence  of  the  Lagrangian  multipliers.  Further  research  into  their  behavior, 
particularly  in  the  presence  of  degeneracy  or  alternate  optimal  solutions,  is  in 
order.  Within  the  context  of  quadratic  programming,  alternate  optimal  solutions 
are  possible  whenever  the  objective  function  has  a  singular  Hessian,  that  is,  the 
Hessian  is  not  positive  definite. 
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APPENDIX  1 


This  appendix  contains  the  proofs  for  the  lemmas  and  theorems  from 
Anstreicher.  which  are  presented  without  proof  in  Chapter  7. 


Lemma  VII.l:  Let  x  be  fe^^sible  in  QP  and  (x,  u)  be  feasible  in  QD.  Then 
f(x)  -d(x.  u)  >  0. 

Proof:  f(x)  —  d{x,  u)  =  b  ^x  +  ~  Since  (x,  u)  is  dual  feasible, 

b^  =  u  —  x^Q, 

f(x)  —  d(x,  u)  =  +  u^Ax  —  x^Qx  +  ^x^Qx  —  u^c 

=  u^(Ax  —  c)  +  jx^Qx  —  x^Qx  +  ^x^Qx 

Now  define  s  =  Ax  —  c. 

f(x)  —  d(x,  u)  =  u^s  +  ~  x^Qx  +  ^x^Qx 
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=  u^s  +  i(x’^  -  x^)Qx +  ^x’^Q(x  -  x) 

=  u'^s  +  |(x  -  X  )'^Qx  -  ^(x  -  X  )’^Qx 
=  ii^s  +  i(x  —  x)  ^Q(x  —  x) 

Thhe  feasibility  of  x  in  QP  implies  s  >  0.  We  also  have  u  >  0  and  Q  positive 
definite.  Thus 

f(x)  -  <hx,  n)  =  u's  +^(x  —  x)^Q(x  -  x)  >  0 


which  proves  the  lemma. 

Lemma  VIL2;  Let  x  6  IR*  and  //  >  0  be  arbitrary.  Then 
i|s~  *PsiP  <||Pi||//  <  S{x.n)-. 

Proof:  llPxllw  =  pI(^Q  +  A'^S  ~ >  pJfA'^S  ~^A)pj.  =  ||s  ~  ’psf-,  for  the  first 
inequality. 

Now  let  u  =  u(x,//)  which  we  define  aa  the  minimizing  value  of  u  resulting  in 
6(x,/x).  Therefore  A^u  =  6  +  Qx,  and 
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=  (2^  -  A  's  -  > )r(lQ  +  A'-S  -  \\)- '(2^  -  A ‘s  -  ‘ ) 

=  - »  ‘  ')r(p2  +  A’'S  -  -A)  -  '(AT(a  -  s  -  ' )) 

.  (lu  _  .)ts  -  lA(iQ  +  A-^S  -  ^'A)-  'A-^S  "  -  e) 

where  a  is  the  maximum  eigenvalue  of 

S-iA(iQ  +  A^S--A)-‘A‘'S-'.  (AM) 

To  evaluate  cr,  we  see  that 

S->A(iQ+ATS-'-A)-'ATS-‘f 

=  S-'A(iQ  + ATS-2A)-'A'"S-‘S-'A(iQ  +  ATS-2A)-‘ATS-‘ 

=  S  -  >  A(iQ  +  A-^S  -  ^a)  -  >(  A^^S  -  2a)(1Q  +  A^S  '  ^a)'  '  A^S  '  '  (A1.2) 

The  maximum  eigenvalue  of  (A1  2)  is  less  than  or  equal  to  that  of  the  following 
matrix  (A1.3) 

S-'A(iQ  +  ATS-^A)-'(J,q  +  ATS-2A)(iQ  +  ATS-2A)-'A'^S''  i' 
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=  S  -  +  A'‘'S  -  '“A)-  1  A^S  -  1 

which  is  (A  1.1).  Since  the  maximum  eigenvalue  of  the  square  of  (A  1.1)  is  less  than 
that  of  (A  1.1),  namely  a.  we  must  have  <t  <  1.  which  implies  ||Pj;||^  <  <5(x./i)“  and 
completes  the  proof  of  the  lemma. 

Lemma  VII. 3:  Given  x  £  !R*,  .suppose  j|S  ~  ^Psjj  <  1.  Let  x  =x-)-p^,  s"  ^s  +  p^, 
and  S  =  diag(s'  ).  Let  H  =  V‘^$(x  .p),  pj  =  —  H  “  ^  V$(x  ,p),  and  =  Ap^  . 
Then  x  €  R*  and 

1  s  ~^Ps  l|<l|Px  III  <^'(X  .p)  <||S“  ^Pslf- <||P^|||  <  (!!'(x,p)- 

Proof:  T  =  s  +  Ps  =  s(e  +  S  ~  ^Ps).  Since  |S  ~  Vs||  <  1>  the  absolute  value  of  each 
element  of  the  vector  S  “  'p^  must  be  less  than  1,  and  each  element  of  (e  +  S  “  ^Ps) 
is  strictly  positive.  The  vector  s  is  also  strictly  positive,  so  s'  must  be  as  well  and 
therefore  x  €  (R*.  Now  consider 

^(x  ,u)  =  min  -nS  u  —  e 

u|  A'"u  =  6  +  Qx  " 

Let  u  satisfy  A^u  =  b  +  Qx.  (One  such  value  is  u  =  u(x,p).)  Then 
X  =  X  —  H  ”  ^  V$(x,p) 

=  x-(lq  +  A^S-2A)-'(9^-AT,-l) 
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=  x-(1q  +  ats--M)->a’-s-'(|!-£) 


Some  manipulation  yields 

lQ(x  -  X- )  =  ^(iQ  +  A^S  -^a)-  'ATS  “  -  e) 

and  using  the  fact  that  +  A^S  ^  =  I  —  A^S  “^A^^Q  +  A^S  "^a)'”  ^ 

vve  obtain 

iQ(x  --  X  )  =  ATS  -  -  ATS  -  2a(1Q  +  aTS  -  ^a)  -  '  aTS  “  '(^  -  e). 

=  ATS-'(^-e)-ATS-''A(lQ  +  ATS-''A)-‘(^-ATs-‘) 

=  ATS  -  '(|i  -  e)  -  ATS  -  ‘(S  -  ‘  A)(iQ  +  ATS  "  ^a)-  -  ATs  “ ') 

=  ATS  -  -  e)  -  ATS  -  '(S  -  ‘a)H  -  ‘  V$(x,p) 

=  ATS-'(|!-e)+ATS-'{s-‘Ap,) 

=  ATS-'(^-e  +  S-‘p.)  (A1.4) 


From  (A1.4)  we  see  that 

Qx  =  Qx  —  A^u  +  pA^S  ~  ’(e  —  S  ”  ^p5) 
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=  —  6  + /iA^S  '(e  —  S  (A1.5) 

so  by  setting  u  =  ~  ^(e  —  S  ~  we  have  A^u  =  6  +  Qx  .  Now  let 

Pg  =  diag(ps)  and  let  (s  ~  ^Ps)^  denote  the  vector  of  the  squares  of  the  elements  of 

S  ~  ^Ps-  Then 

6(x  ,n)  =  min  iSu  — e 

u|  A'^u^ft  +  Qx  ' 

<  j  u  —  e  I 

=  (S  +  PJS~’(e-S"‘ps)-e 

=  ||(l  +  S-‘pJe-S-‘p,)-e| 

=  le-S"‘pj  +  S''ps-S“'PsS"'ps-e]| 

=  |l(s-‘pd11 
<|s-‘p»IP 

This  proves  S(x  ,fx)  <  ||s  “  ^PsP»  and  the  remaining  inequalities  are  a  result  of 
Lemma  VII. 2. 

2 

Lemma  VIL4:  If  x  G  (R*  and  ||Px|(ff  <  then  $(x,//)  —  ^(x(f2).fi)  <  ■■ 
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Proof:  For  an  arbitrary  fixed  /z  >  0,  $(  •  ,fi)  is  convex,  so 

$(x,/z)  -  $(x  +  Pj„/i)  <  -  V$(x,/i)'^Pj.  =  ||P^||5^  (A1.6) 

Let  X®  =  X,  and  define  |x^,  x^.  x“,  ...  |  as  the  sequence  of  points  obtained  by 
generating  a  succession  of  NMvton  steps  starting  at  x.  From  (A1.6)  and  Lemma 
VIL3. 

oo 

$(x,/z)  -  <I'(x(/i),/i)  =  ^($(x',^)  -  $(x*  +  \/z)) 

z=0 


Vf(x)'^p^  <  f(x  +  p^)  -  f(x)  <  Vf(x  +  p^)'^Pi. 


(A1.7) 
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If  the  right  hand  and  left  hand  sides  of  (Al.7)  can  be  appropriately  bounded,  the 
lemma  can  be  proved.  Beginning  with  the  left  hand  side  and  the  fact  that 
Vf(x)  =  Qx  +  b  =  pV$(x,/i)  +  pA^S  “  we  see 

Vf(x)'^p^  =  pV$(x.p)'^p^  +  He'^S  ~  ^Apj. 

=  /xV$(x,/i)''pj.  +  /ie'^S“^Ps.  (A1.8) 

Xow  V$(x,p)'^p^  =  V^(x4ip'{  -  H  -  ^V$(x.//))  =  -  V$(x,//)TH  -  “  ‘V^>(x,/i) 

=  —11  Pi  II//-  From  the  Cauchy-Schwartz  Inequality  and  Lemma  VII. 2  we  have 

e'^S~^Ps>  -|e'^S"’p5|>  -||e||||S“  Vs||>  - 'Im  ||Pi||//- 

Substitution  into  (A1.8)  yields 

Vf(x)'^p^>  -/i||p^l|^-//>fm|lp^|j;^.  (A1.9) 

Using  (A1.5)  and  noting  that  x  =  x  +  pj.,  we  can  bound  the  right  hand  side  of 
(Al.7)  as  follows: 

Vf(x  +  Pj)'^Pj  =  (  Qx  +  6)'^Pj. 

=  /xpjA'^S"*(e-S~^Ps) 
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=  /ie'^(s  “  ^A.pj.)  -  /ipJ'A'^S  ~  ’S  “  Vs 

= /xe'^(s  “  Vs)  -  J  S  “  Vs  IP 

</ZAim||Px||yy.  (ALIO) 


We  can  now  subsitute  (A1.9)  and  (Al.lO)  into  (A1.8)  to  achieve 


~  mxlH  -  llPxll//  <  f(x  +  Pi) -f(x)  < /iNfm  1|Pj.||^ 

-  /iVm  II  Pj.  ||5/  -  /x^^m  II  p^  11;^  <  f(x  +  pj.)  -  f(x)  <  ||  p^  ||^  +  /i^  ||  p^ 


lf(x  +  Px)  -  f(x)  I  <  ||Pilli/(l  +|lPx||//)  (ALII) 

Again  let  x®  =  x  and  {x®,  x^  x^,  .  .  .}  denote  the  sequence  of  points  obtained  by 
generating  a  succession  of  Newton  steps  starting  at  x.  Applying  (ALII)  and 
Lemma  VII. 3  we  have 


|f(x)-f(x(M))l  = 


£(f(x')-f(x‘  +  ')) 

t=0 


<  £|f(x‘)-f(x'  +  ')| 
1=0 
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Since  (1  + 11  Pi  IIt/  )  ^  1  + 1|  Pi  \\~u  J  for  all  i  >  0  this  becomes 


oo  j 

|f(x)  -f(x(/^))|<  /i^^rn(l  +||Pi||^)X!llPx||// 

i=0 


-  l-llPxllff 


which  proves  the  lemma. 

Theorem  VII.l:  Let  z*  denote  the  minimum  value  of  f(x)  for  problem  QP.  After 
0(L/^)  outer  iterations,  the  algorithm  obtains  such  that  mfi^  <  2  ~  an 

a  point  X  =  x^  *  such  that  f(x)  —  z*  <  ^1  +  2  ~ 

Proof:  The  duality  gap  at  x(/i^)  is  and  at  the  end  of  the  outer  iteration 

we  have  |(Pi||/^  <  t.  This  along  with  Lemma  VII. 5  gives  us 

f(x)  -z*  =  f(x(/))  -z*  +  f(x)  -  f(x{/t^))  <  m/i*^  + 

S  "■»*'{ l+ffel^)  (A1.I5 

To  show  mfi^  <  2  note  that  m/x^  =  m(l  —  and  <  2^^^^  is  assumed. 

So  we  need  a  value  of  k  such  that 


Taking  logarithms, 

kMl-^)  +  0(L)+/n(m)<  -0(L) 

,  >  Q(L)  +  /n(m) 

-  -Hi-e) 

and  since  ln[m)  <  0(L)  and  6  <  — /n(l  —  ^),  we  obi  i 

k>0(|) 


which  completes  the  proof. 


Lemma  VII.6: 


Let  A  = 


1 


1 


Then 


$(x,^)  -${x  +  A  p^,/i)  >||Pjp||y^-/n(l  +||Px||^) 


Proof:  The  Taylor  series  expansion  for  $(x  +  Ap^.,;/)  is 

$(x  4-  Ap^,/^)  =  $(x,/z)  +  AV$(x,^)'^p^  +  ^A^pjHpj.  +  Yl  t^-  (A1.13) 

j=3 

where  tj  denotes  the  term  of  the  Taylor  expansion.  For  j  >  3,  the  term  may 
be  computed  as 
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2  =  1 


where  (s  ^Ps)f  denotes  the  power  of  the  component  of  S  ^p,.  Thus 


|f,|<^g|(s-'pJ|<^(E|(s-'p.), 

■'2  =  1  \2=1 


W^  =  ^|s-‘P.IF 


and  bv  Lemma  VII. 2, 


J 


(A1.14) 


Expanding  the  linear  and  quadratic  terms  of  the  Taylor  series  and  again  using 
Lemma  VIL2.  we  obtain 

AV$(x,/i)'^p^  +  iA2pjHp^  =  -  XW^x,^ifE  -  iRH  "  ^  V$(x,p)  +  iA2pjHp, 


=  (_A  +  iA2)||p^ll^. 


(A1.15) 


Substituting  (A1.14)  and  (A1.15)  into  (Al.13)  yields 


oo 


(A1.16) 


$(x  + Ap^,/2)  < 

J=^ 

Furthermore,  the  Taylor  series  expansion  of  —  ln(^l  —  A||pj.||^  j  is 


-  ^1  ~  A||Pj.||^]  -  -  Ml)  + -^llPill// +  +  l{'^llPxll^)‘^  +  .  .  . 
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which  implies 


oo  >  j  . 

^yl|Px|li/=  “  ~  ^llPxlltf  )- -^llPxll^  -  2{'^!|Pi|li/)^  (A1.17) 

Now  we  substitute  (A1.17)  into  (Al.16)  and  find 

$(x  +  Ap^,^)  <  -  A||Pj.||j^  -  /n(l  -  A||Pj.||^)-  A||p^ll;^. 

The  right  hand  side  is  minimized  bv  A  =  A  =  — 5^ - ,  so 

1+llPxll/f’ 

$(x,/i)  -$(x  +  A  Px,/i)  >  A||p^||^(l  +||p^||^^)-ln(l  -AIIP^II;/) 

=  llPx||ff-^l+||Px||//) 

which  proves  the  lemma. 

Corollary  VII.6:  Let  r  <  l|Pi||//  <  Then  a  pure  Newton  step  (A  =  1)  obtains 

$(x,/x)  —  $(x  +  PxiH)  >  A  =  +  ln{l  —  r)  +  r  >  0 

Proof:  Using  a  step  size  of  A  =  1,  from  Lemma  VIL6  we  have 
$(x  +  p^,/z)  <${x,/x)  -||p^||^_/n(l-||p^||^)_||p_^||^ 
or  equivalently 
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-$(x  +  p^,/y)  >||p^||^  +  /n(l  -|jPr||//)  +  ||Px||//- 


(A1.18) 


The  right  hand  side  of  (A1.18)  is  positive  and  increasing  for  0  <  jlPxlly/  <  so 
0  <  r  <  II  Pj.  11^^  <  i  implies 

0  <  A  =  T“  +  /n(l  —  r)  +  r  <  $(x,p)  -  $(x  +  Px,/^) 


which  completes  the  corollary. 


r(l  +  r) 


Theorem  VII.2:  Let  A  =  r  —  ln{l  +  r),  E*  =  — - — ij,  and  N  denote 

the  number  of  inner  iterations  required  for  an  arbitrary  outer  iteratation.  Then  N 


is  bounded  as: 


N  <^j-^(^m  +  2Ey>|m)  +  E<^^ 


Proof:  Anstreicher’s  proof  is  a  generalization  of  a  proof  by  Gonzaga  for  the  linear 
case.  The  (k  +  1)®^  outer  iteration  begins  with  the  point  (x^.p^)  and 
|p(x^.//^“  <  7-  for  H  =  H(x*^,/i^~  ').  Let  N  denote  the  number  of  inner 

iterations.  Lemma  VII. 6  shows  that  each  inner  iteration  will  decrease  $(  •  ,//^)  by 
at  least  A.  Following  the  N  iterations,  we  have  x^’  ^  with  ||p(x^'*'  ^P^)x|jy  ^ 

H  =  H(x^  +  ^/).  So 


AN  <  $(x  V)  -  $(x*^  +  ^^*). 


(A1.19) 
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The  right  hand  side  can  be  bounded.  From  the  definition  of 


and  this  leads  to 


$(x^/)  -  ‘h(x^  +  =  $(x^/-  1)  -  $(x^  +  ~  M 


Analyzing  (f(x^’)  —  f(x^  we  use  the  the  fact  that  and  x^  ^  are 
approximately  centered  with  respect  to  /i*“  ^  and  respectively,  and  Lemma 
VII. 5  to  obtain 


f(x*)  -  +  ')  <  f(x(,i‘  -  '))  +  -  'Vm  -  f(x(/.‘))  + 

=  f(x(^‘  -  > ))  -  f(x(/))  +  E/(2  -  «)//*  -  '.im. 

Recall  the  dual  objective  function  d(x,u)  from  Lemma  VII.  1,  and  u(x,yi)  from 
Lemma  VII. 2.  Since  u(x(/z),/r)  is  dual  fea.sible,  with  the  monotonicity  with  respect 
to  fi  of  the  dual  objective  function  we  have 

f(x*)  -  f(x‘  +  >)  <  f(x(/<‘  -  ‘ -  d(x(/  -  ‘ u(x(„‘  -  ' ) 
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+  d(x(//-).  U(x(/)./))  -  f(x(/i^-))  +  E^(2  -  «)/  - 


<f(x(/-h)-(l(x(/-M-  u(x(/-M./  ^)) 

-  (f(x(/))  -  a(x(/‘).  u(x(/)./’)))  +  2EyM^'  - 

<  +  2E ~  ^^Jm 

<  "  '  +  2Ey^^'  "  W 

</-^(^m  +  2E^>[m|  (A1.21) 

We  can  now  use  Lemnaa  VII. 4  with  ||Pxll//  ^  nearly  centered  with  respect  to 
*,  and  x(/'  ')  minimizing  to  achieve 

<  $(x*,/'  -$(x(/“ 

<E.  (AI.22) 


Finally,  substitution  of  (A1.21)  and  (A1.22)  into  (Al.20)  yields 
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-  $(x^'  +  ‘  ./■)  <  +  2E^^) 

<  [Om  +  2Ejm)  +  E^ 


and  substitution  of  this  inequality  into  (A1.19)  proves  the  theorem. 
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